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Abstract. The three-body problem, which describes three masses interacting 
through Newtonian gravity without any restrictions imposed on the initial positions 
and velocities of these masses, has attracted the attention of many scientists for more 
than 300 years. In this paper, we present a review of the three-body problem in 
the context of both historical and modern developments. We describe the general 
and restricted (circular and elliptic) three-body problems, different analytical and 
numerical methods of finding solutions, methods for performing stability analysis, 
search for periodic orbits and resonances, and application of the results to some 
interesting astronomical and space dynamical settings. We also provide a brief 
presentation of the general and restricted relativistic three-body problems, and discuss 
their astronomical applications. 
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In the three-body problem, three bodies move in space under their mutual gravitational 
interactions as described by Newton’s theory of gravity. Solutions of this problem require 
that future and past motions of the bodies be uniquely determined based solely on their 
present positions and velocities. In general, the motions of the bodies take place in 
three dimensions (3D), and there are no restrictions on their masses nor on the initial 
conditions. Thus, we refer to this as the general three-body problem. At hrst glance, the 
difficulty of the problem is not obvious, especially when considering that the two-body 
problem has well-known closed form solutions given in terms of elementary functions. 
Adding one extra body makes the problem too complicated to obtain similar types of 
solutions. In the past, many physicists, astronomers and mathematicians attempted 
unsuccessfully to hnd closed form solutions to the three-body problem. Such solutions 
do not exist because motions of the three bodies are in general unpredictable, which 
makes the three-body problem one of the most challenging problems in the history of 
science. 

In celestial mechanics, the general three-body problem deals with gravitationally 
interacting astronomical bodies and intends to predict their motions. In our Solar 
System, the planets and asteroids move around the Sun, while the moons orbit their 
host planets, which in turn also move around the Sun. As typical examples of the 
three-body problem, we may consider the Sun-planet-planet, Sun-planet-moon, or Sun- 
planet-asteroid systems. The three-body problem representing the latter system can 
be signihcantly simplihed because the mass of the asteroid is always negligible when 
compared to the mass of either the Sun or the planet, which means that the gravitational 
influence of the asteroid on the planet and the Sun can be omitted from the theory. If 
this condition is satished, then the general three-body problem becomes the restricted 
three-body problem, and there are two possibilities, namely, the two bodies with dominant 
masses move around their center of mass either along circular or elliptic orbits, which 
leads to the respective circular or elliptic restricted three-body problems. It is the circular 
restricted three-body problem (CR3BP) that has been the most extensively studied. 

The three-body problem has been studied for over three hundred years. We wish to 
provide a brief historical overview of the most signihcant results that have contributed 
to the current developments in the held. We divide this overview into two parts: hrst, 
we highlight the time period from the publication of Newton’s Principia in 1687 to the 
publication of Poincare’s Les Methodes Nouvelles de la Mecanique Celeste published 
in 1892 (volume 1 and 2) and in 1899 (volume 3); in the second part we focus on 
developments that took place from 1900 to the present time. 


1.1. From Newton to Poincare 


The three-body problem was formulated and studied by lNewtonl 1687l | in his Principia, 
where he considered the motion of the Earth and the Moon around the Sun. In this 
problem there are key elements such as the ratio of the Moon’s mass. Mm, to the 
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Earth’s mass, Me, that is Mm/Me = 0.0123, which is a small but not negligible 
number. Another element is the tilt of the Moon’s orbit to the Earth’s orbit, about 
5°. This makes the problem completely predictable despite the fact that the Earth’s 
orbit around the Sun and the Moon’s orbit around the Earth are almost perfect circles. 
Newton stated that the problem was very difficult to solve; however, he was able to 
obtain an approximate solution which agreed with observations to within 8%. 

A special form of the general three-body problem was proposed by lEulerl 1767 . 


He considered three bodies of arbitrary (hnite) masses and placed them along a straight 
line. Euler showed that the bodies would always stay on this line for suitable initial 
conditions, and that the line would rotate about the center of mass of the bodi es, leading 
to per iodic motions of all three bodies along ellipses. Around the same time, iLagrange 
1772l | found a second class of periodic orbits in the general three-body problem. He 
showed that if the bodies were positioned in such a way that they form a triangle of 
equal sides which would move along ellipses for certain initial conditions, preserving 
always their original conhguration. The Euler and Lagrange solutions are now known 
as particular solutions to the general three-body problem, and they will be further 
disc ussed in Sect ion 3. 

Euler! 17671 was the hrst to formulate the CR3BP in a rotating (or synodic) 
coordinate system. This was an important development in the study of the three- 
body problem. Lagrange also studied the CR3BP and demonstrated that there were 
hve equilibrium points (now known as the Lagrange points) at which the gravitational 
forces of the bodies balanced out; the Trojan asteroids discovered in 1906 along Jupiter’s 
orbit occupy space close to two of the Lagrange points. 

who used 
to 


836 


Euler 1767 


Important contributions to the CR3BP were made by I Jacobi 
the synodic (rotating) coordinate system originally introduced by 
demonstrate that there was an in t egral of mo tion, which is now named after him. The 
Jacobi integral was used by Hill 1877 . Il878| to determine the motion of an asteroid 
in the three-body problem and to introduce the so-called zero velocity curves (ZVC), 
which establish regions in space where the bodies are allowed to move. Hill considered a 
special case of the CR3BP in which two masses were much smaller than the hrst one (the 
problem is now known as the Hill problem), and in this way he discovered a new class 
of periodic solutions. His main contribution was to present a new approach to solve the 
Sun-Earth-Moon three-body pr oblem. After a 


formulation of the problem by 

Newton 

1687 

, Hill developed his lunar theory, which 

with some modihcations made bv 

Brown! 

1896 

, is still being used today in celestial 


mechanics Gutzwillerl . Il998 


In the second half of the nineteenth century, Poincare studied and advanced the 
solution of the three-body problem. His monumental three-volume book Les Methodes 
Nouvelles de la Mecanique Celeste, originally published in 1892-99, contains his most 
important contribution to the study of the CR3BP; the book was translated as New 
Methods of Celestial Mechanics and edited by L. D. Goroff in 1993. In this work, 
Poincare developed a number of new qualitative methods to solve differential equations 
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and used them to identify and study possible periodic orbits, while demonstrating non- 
integrability of the system of equations describing the three-body problem. Poincare’s 
new methods allowed him to identify the unpredictability of the problem, and to discover 
the hrst manifestation of a new phenomenon, which is now commonly known as chaos. 
Poincare submitted some of those results for the King Oscar II Birthday Competition 
and he was a warded the prize. Th e competition and Poincar e’s ot her contr i bution s are 

and 


described by iBarrow-GreenI 1997 
will not be repeated here. 


Diacu and Holmes 


1996 


and Peterson 


1993 


1.2. From Poincare to the present time 


Poincar 

e 

1892 

studied 

orbits ! 

ad, 

1877, 

1878 


s tudied Hill’s problem and generalized Hill’s dehnition of periodic 
He was able to choose such initial conditions that resulted 
in periodic orbits in the special restricted three-body problem. Poincare’s work on 
the existence of periodic so l utions in dynamical systems with one degree of freedom 


f 


was extended by iBendixsonI [1901l | . He formulated and proved a theorem that is now 


known as the PoincarGBendixson theorem, which gives a criterion for the existence of 
periodic solutions in such systems. The work of Poincare led to syste matic searches 
for p e riodic orbits in t he th ree-b ody problem and their classihca tion by iDarwinI 18971 . 


19091 , iMoulton et al.l 192nl | and IStromgren and Pedersen 


stu dies of stabi li ty of such orbits, which was initiated bv iPoincara |1892| then continued 


19221. as well as extensive 


by Levi-Civita 1901] and Lvapnnov jl9n7j |. who signihcantly generalized Poincare 


approach and his results. 

Further generalizations a nd extensions of Poincare’s ideas on the stability of 


periodic orbits were pursued bv iBirkhoffl 1912| . He introduced the concept of recurrent 


motion, which requires a sufficiently long time so that the motion comes arbitrarily 
close to all its states o f motion, and showed its relationship to orbital stability. In 
addition, iBirkhoffl 1913| proved the ‘Last Geometric Theorem’ formulated by Poincare. 


The Poincare-Birkhoff theorem states that there are inhnitely many periodic or bits near 
any st able periodic orbit; it also implies the existence of quasi-periodic orbits. iBirkhoff 
1915l | also developed a topological model for the restricted three-body problem in which 
an asteroid was conhned to move inside an oval about one of the larger masses. 

By searching for periodic and non-periodic solutions to the three-body problem, 
researchers realized that the differential equations describing the problem contain 
singularities at which the solutions are abruptly terminated. An example of such a 
termination are collisions between two or even all three bodies; thus, we have the so- 
called collision singularities. Once their presence is established, the next challenging 
task is to eliminate them by a process called regulariz ation. A n imp o rtant work 
on singularities in the three-body problem was done by Painleve 1896 . Il897j |. He 


determined that collisions were the only singularities and that collisions can be excluded 
by setting certain initial conditions for which the equations of motion of the three-body 
problem could be integrable using power series solutions. Painleve did not hnd such 
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solutions but his work had a strong impact on other researchers in the held. The 
pr oblem of the regula r ization of collisio ns bet ween two b odies was investig ated mainly 


by Ihevi-Civital 19031 , iBisconcinil lOnoj ] , and ISnndmanI 19071 . Il909l Il912j | . They also 


studied the triple c ollisions and formulated theorems that allowed establishing conditions 


for such collisions Siegel and Moserl . Il991 


Realizing that the possibility of hnding closed-form solutions to the t hree-body 
proble m was unlikely, researchers considered hnding inhnite series solutions. iDelannav 
1867| used canonical varia bles in the pertu rbat ion theory of the three-body problem. 
Following Delaunay’s work. Ihindstedtl |l884| and lGvldenl jl893| introduced inhnite series 
and investigated the rates of the ir convergence, but they were not able to present a 
formal proof of the convergence. iPoincara 1892| clarihed the concepts of convergence 


and divergence of the inhnit e series and exposed certain haws in the methods developed 
by these authors. IPainleva |l89fil . Il897l | also attempted to hnd solutions to the three- 
body problem given in terms of inhnite series but failed as the others did before him. 
However, he clearly stated that such s olutions w e re pos sible, at least in principle. 


Indeed, Painleve was right and ISnndman 


1912 


found a complete solution to 
the three-body problem given in terms of a power series expansion. Unfortunately, 
Sandman’s solution converges very slowly so that it cannot be used for any practical 
applications. An important question to ask is: Does the existence of Sandman’s solution 
contradic t the unpredictability of motions in the three-body problem postulated by 


Poincare 


1892 


? Well, the answer is that it does not because trajectories of any of 
the three bodies cannot be determined directly from Sandman’s solution, which means 
that the trajectories can still be fully unpredictable, as this is observed in numerical 
simulations (see Section 4, 5 and 8). 



integrable and non-integrable. As shown above, even the systems that are in principle 
non-integrable may still have periodic solutions (orbits) depending on a set of initial 
conditions. Poincare, Birkhoff and others considered quasi-periodic solutions (orbits) in 
such systems. The very fundamental question they tried to answer was: What happens 
to solutions of an integrable Hamiltonian system when the governing e quations are 
slight ly perturbed? The correct answer to this question was hrst given by iKolmogorov 


1954 1 but withou t a fo rmal proof, which was later supplied independently by 


Moser 


1962 1 and Arnold 1963| ; the formal theorem is now known as the Kolmogorov-Arnold- 


Moser (KAM) theorem. This theorem plays an important role in the three-body problem 
and in other Hamiltonian dynamical systems Hilbornl. Il994 . 


m 


As mentioned above, periodic orbits in the three-body problem were discovered 

By utilizing computer simulations. 


Henon 1965. 

1974 

and 

Szebehelv 

1967 


1967(1 found many p eriodic orbits and classihed 


them. Recent work by ISuvakov and Dmitrasinovid 2013| shows that periodic orbits 
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are still being discovered. An interesting problem was investigated numerically 


by ISzebehelv and PetersI 1967alJb| . who considered three objects with their masses 


proportional to 3, 4, and 5 located at the vertices of a Pythagorean triangle with 
the sides equal to 3, 4, and 5 length units; they were able to show that two of 
these objects form a binary system but the third one was expelled at high speed 


Moore. 

1993, 

Chenciner and Monteomerv. 

2000, 

Suvakov and Dmitrasinovic. 

2013 


and 

13 new periodic solutions for th e general three-body problem with equal masses 
Suvakov and Dmitrasinovig . 12013 . 


Modern applications of the three-body problem have been greatly extended to 
include the Earth, Moon, and artihcial satellites, as well as the recently discovered 
distant extrasolar planetary systems (exoplanets). With the recent progress in detection 
techniques and instrumentation, there are over a thousand confirmed exoplanets orbiting 
single stars, either both or a single component within binary stars, and even triple 
stellar system^ along with thousands of exoplanet candidates identihed by the NASA’s 
Kepler space telescop^. A typical three-body system would be a single star hosting two 
exoplanets or a binary stellar system hosting one exoplanet; we discuss such systems in 
this paper. 


Whittaker 

19371, 

Winl 

tnerl 

194ll, 

Pollard 19f 

iol. Dan 

bvl 198? 

^ . Siesrel and Moser 

1991 . 

Murrav and DermottI 

1999 

. Rov 

2 OO 5 I 1 and 

Beutler 

20G 

51, 

and in review papers such 

as Holmes 

1990 . 

Szebeheb 

/ 199? 


diitzwille 

rl 19981 

and 

Ito and Tanikawa 

2007 

|. The 

books bv 

r’oinca] 

4 1892 . 

Szebeheb 

[iMZl, 

Marcha 

19901 

Barrow-Green 

199 

7 and 


Valtonen and KarttunenI 2006j | are devoted exclusively to the three-body problem and 
also discuss its many applications. 

This paper is organized as follows: in Section 2 we formulate and discuss the general 
three-body problem; then we describe analytical and numerical studies of this problem 
in Sections 3 and 4, respectively; applications of the general three-body problem are 
presented and discussed in Section 5; the circular and elliptic restricted three-body are 
formulated and discussed in Sections 6 and 7, respectively; applications of the restricted 
three-body problem are given in Section 8; the relativistic three-body problem and its 
astrophysical applications are described in Section 9; and we end with a summary and 
concluding remarks in Section 10. 


2. The general three-body problem 

2.1. Basic formulation 

In the general three-body problem, three bodies of arbitrary masses move in a three 
dimensional (3D) space under their mutual gravitational interactions; however, if the 

f http://exoplanet.eu/catalog/ 

§ http://kepler.nasa.gov/mission/ciiscoveries/candidates/ 
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motions of the bodies are confined to one plane the problem is called the planar general 
three-body problem. Throughout this paper, we use Newton’s theory of gravity to 
describe the gravitational interactions between the three bodies; the only exception 
is Section 9 in which we discuss the relativistic three-body problem. 

To fully solve the general three-body problem, it is required for the past and future 
motions of the bodies to be uniquely determined by their present positions and velocities. 
Let us denote the three masses by Mj, where i = 1, 2 and 3, and their positions with 
respect to the origin of an inertial Cartesian coordinate system by the vectors Ri, and 
dehne the position of one body with respect to another by = Rj—Ri, where Tij = —rji, 
j = 1, 2, 3 and i ^ j. With Newton’s gravitational force being the only force acting 
upon the bodies, the resulting equations of motion are 


M, 


d^Ri 

dR 


gE 


MM,- 


IJ 5 


( 1 ) 


where G is the universal gravitational constant. 

The above set of 3 mutually coupled, second-order, ordinary differential equations 
(ODEs) can be written explicitly in terms of the components of the vec tor Ri, 
which means that there are 9 second-order ODEs; see iBroucke and LassI 1973| for an 
interesting and useful form of the equations. Since these equations can be converted into 
sets of 2 hrst-order ODEs, there are actually 18 hrst-order ODEs to fully describe the 
general three-body problem. A standard mathematical procedure Whittaker! . Il937 
allows solving such ODEs by quadratures if independent integrals of motion exist. 
Therefore, we shall now determine the n umber of integrals of motion for the general 
three-body problem jMcCord et al.l . Il998 . 


2.2. Integrals of motion 


Let us sum Eq. ([T]) over all three bodies and take into account the symmetry condition 
Tij = —rji. The result is 


3 


Y,M. 


d'^Rj 

dR 


0 , 


( 2 ) 


and after integration, we obtain 



2 = 1 


(3) 


where Ci = const. One more integration yields 

3 

MiR, = CR + C 2 , (4) 

2=1 


with C 2 = const. 


















The three-body problem 


^ 3^3 

Since the center of mass is defined as Rem = X] (0!) determines 

i=l i=l 

its motion, while Eq. (Ej) shows that it moves with a constant velocity. The vectors Ci 
and C 2 are the integrals of motion, thus we have 6 integrals of motion by taking the 
components of these vectors. 

The conservation of angular momentum around the center of the inertial Cartesian 
coordinate system in the general three-body problem gives other integrals of motion. 
To show this, we take a vector product of Ri with Eq. ([2]), and obtain 


<PRi 

^ MiRi X = 0 , 

i=l 


dR 


which after integration yields 


x-^ = C'3, 


2=1 


dt 


(5) 


( 6 ) 


where C 3 = const. Hence, there are 3 more integrals of motion. 

An additional integral of motion is related to the conservation of the total energy 
of the system. With the kinetic energy given by 



2=1 


dRi 

dt 


(7) 


and the potential energy Ep^t defined as 


Epot - 2^ 2^- 


i=l j=l 


r. 




( 8 ) 


where i j, we have the total energy Etot = Ei^in + Epot = =const to be also an 
integral of motion. 

According to the above results, there are 10 (classical) integrals of motion in the 
general three-body problem, which means that the set of 18 first-order ODEs can be 
reduced to 8 first-order ODEs. Actually, there are 2 other integrals of motion, one 
related to the elimination of time and the other to the elimination of the so-called 
ascending node. The time can be eli minated by tran s forming one of the d ep endent 


variab les as an independent variable Szebehelyi . 119671 . IB arrow- Greenl . Il997l |. I.Tacobi 


1843| changed variables so that two bodies would orbit a third one, and showed that the 


difference in longitud e between the a s cendi ng nodes is fixed at tt radians, which becomes 
an integral of motion Barrow-Green . 1997j |. Having obtained the 12 integrals of motion. 


other independent integrals of motion exist Brims, 

1887, 

Poincare. 

1892, 

Whittaker. 

1937, 

Szebehelv. 

1967, 

Valtonen and Karttunen. 

2006 
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Since i^kin > 0 and Epot < 0 (see Eqs ([7]) ([HD, respectively), i?tot or C 4 can either be 

positive, negat ive o r zero, and 6*4 can be used to classify motions of the general three- 

body problem jRovl . I2nn5l |. In the case C 4 > 0, the three-body system must split, which 

means that one body is ejected while the remaining two bodies form a binary system. 

The special case of (^4 = 0 is unlikely to take place in Nature; however, if it indeed 

occurred, it would result in one body escaping the system. Finally, the case C 4 < 0 may 

lead to either escape or periodic orbits with the result depending on the value of the 

3 


moment of inertia given by J = ^ for details see iRov 

Valtonen and Karttunen bnnfil .' 


2005 


Marchal 1990 and 


In addition to the integrals of motions, the virial theorem < Fkin > = - < Epot > /2, 
where < Fkin > and < Epot > is the time average kinetic and potential energy, 
respectively, can be used to determine the stability of the three-body problem and its 


statistical properties Valtonen and Karttunen! . l2006| . The system is unstable if its time 


average kinetic energy is more than two times higher than its time average potential 
energy. 


2.3. Another formulation 

The standard formulation of the general three-body problem presented in the previous 
section is often replaced by either the Hamiltonian formulation extensively used by 
Poincare 


1892 


or by the variational formulation described in detail and used by 


Siegel and Moser! 199l| . Here, we present only the Hamiltonian formulation. 


Let us consider the so-called natural units and introduce G = 1 in Eqs ([ID and ([HD. 
Moreover, we write Ri = {Ru, R 2 i, Rsi) = Qki, where Ru, R 2 i and R^i are components 
of the vector R^ in the inertial Cartesian coordinate system, and fc = 1, 2 and 3. Using 
this notation, we dehne the momentum, pki, as 


Pki 


dqki 

dt 


and the kinetic energy as 


Eto. = E """ 


k,i=l 


2M, ’ 


(9) 


( 10 ) 


and introduce the Hamiltonian H = Ekin + Epot that allows us to write the equations 
of motion in the following Hamiltonian form 


dqki _ dH 
dt dpki 


and 


dpki 

dt 


dH 

dqki 


( 11 ) 


which is again the set of 18 hrst-order ODEs equivalent to the set of 18 hrst-order ODEs 
given by Eq. o. 

The equations of motion describing the general three-body problem (either Eq. ([ID 
or Eq. ([HD) can be further reduced by considering the general Hill problem, in which 
Ml » M 2 and Mi » M 3 , however, M 3 is not negligible when compared to M 2 . 
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The governing equation s describing this problem are derived by ISzebehelvi 1967l | and 


Simp and Stuchi 


2000 


who also discussed their applications to the Sun-Earth-Moon 
problem, known as the Hill lunar theory. In this paper, we limit our discussion of the 
Hill problem to its circular and elliptic versions of the restricted three-body problem, 
and describe them in Section 6.5 and 7.3, respectively. 


2 . 4 . General properties of solutions of ODEs 

Solutions of the ODEs describing the general three-body problem can be represented 
by curves in 3D. In a special case, a solution can be a single point known also as a 
fixed point or as an equilibrium solution; a trajectory can either reach the fixed point 
or approach it asymptotically. Typically, periodic orbits are centered around the hxed 
points, which are called stable points or stable centers. Moreover, if a trajectory spirals 
toward a hxed point or moves away from it, the point is called a spiral sink or spiral 
source, respectively. There can also be a saddle point in which two trajectories approach 
the point and two leave it, but all other trajectories are kept away from it. 

From a mathematical point of view, it is required that solutions to ODEs exist and 
that they are unique. The existence can be either global when a solution is dehned for 
any time in the past, present and future, or local when it is dehned only for a short period 
of time. The uniqueness of solution means that there is only one solution at each point. 
For given ODEs, the existence and uniqueness are typically stated by mathematical 
theorems. Let us consider the following hrst-order ODE y'{x) = f{y{x),x), where 
y' = dy/dx, and the intial condition is y{xo) = j/o, with xq being the initial value of 
xe[xQ — e,xo e\. Picard’s (or Lipschitz and Cauchy’s) existence theorem states that if 
f{y{x),x) is a Lipschitz continuous function in y and x, then there is £ > 0 such that a 
unique solution y{x) exists on the interval [xq — e,XQ + e\. 

Another requirement is that the problem is well-dehned, which means that a 
solution must be continuous with respect to the initial data. Being aware of the above 
mathematical requirements, Poincare studied the problem of solving ODEs from another 
perspective by developing q ualitative methods, which he described in the third volume 
of his book Poincare . Il892l |: Poincare’s work was done mainly for the CR3BP, and we 
describe it in Section 6. 


3. Analytical studies of the general three-body problem 

3.1. Euler and Lagrange periodic solutions 

As already mentioned in Section 1, two differ ent classes o f per iodic solut i ons t o the 
general three-body problem were obtained by lEulerl 1767l | and iLagrangd 1772l | who 
considered a linear and triangle conhguration of the three bodies, respectively, and 
demonstrated that in both cases the bodies move along elliptic orbits. Since the masses 
of the three bodies are arbitrary (finite), the problem is considered here as the general 
three-body problem despite the limitation on the bodies mutual positions; because of 
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Figure 1. Illustration of the Euler solution with a line joining three masses Mi, M 2 
and M 3 located at the corresponding points A, B, and C. 


this limitation, the Euler and Lagrange solutions are called the particular solutions. 
The existence of the special conhgurations (known as central conhgurations in celestial 
mechanics) plays an important role in studies of orbital stabil ity. Detailed derivations 


of the Euler and Lagrange solutions are given by iDanbv 


1988 


In the Euler solution, the initial linear conhguration of the three bodies is 
maintained (it becomes the central conhguration), if the ratio AB/BC (see Fig. [1]) 
has a certain value that depends on the masses, and if suitable initial conditions are 
specihed. Euler proved that the line AC would rotate about the center of mass of 
the bodies leading to periodic motions of all three bodies along ellipses. Moreover, he 
also demonstrated that the ratio AB/BC would remain unchanged along AC during the 
motion of the bodies. Since the three bodies can be ordered in three different ways along 
the line, there are three solutions corresponding to the ordering of the bodies. However, 
it must be noted that the Euler solution is unstable against small displacements. 

Now, in the Lagrange solution, the initial conhguration is an equilateral triangle 
and the three bodies are located at its vertices. Lagrange proved that for suitable initial 
conditions, the initial conhguration is maintained (becomes the central conhguration) 
and that the orbits of the three bodies remain elliptical for the duration of the motion 
(see Fig. ED. Despite the fact that the central conhguration is preserved, and the 
triangle changes its size and orientation as the bodies move, the triangle remains 
equilateral. Since the triangle can be oriented in two diherent ways, there are two 
solutions corresponding to the orientation of the tri angle. Regions of stability and 


instability of the Lagrange solutions were identihed by iMansillal 2006 


3.2. Other periodic solutions 


Hill] 1877L Il878| found periodic o r bits i n the problem that is now known as the Hill 
problem (see Section 1). Poincare jl892 | developed a method that allowed him to hnd 
periodic orbits in the CR3BP (see Section 6). However, the original version of his 
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c 



Figure 2. Illustration of the Lagrange solution with an equilateral triangle joining 
three masses Mi, M 2 and M 3 located at the corresponding points A, B, and C. 


metho d can also be applied to the general three-body problem. We follow iPoincare 


1892| and consider the autonomous Hamiltonian system given by 


dxi 

dt 


dF 

dVi 


and 


dyi 

dt 


OF 

dxi ’ 


( 12 ) 


and F = Fq + pFi -|- p,F 2 where Fq = Fq{x), however, Fi, F 2 , ... are functions of 
both X and y and they are periodic (with a period of 27r) with respect to y. Moreover, 
/i is a parameter that depends on mass, where z = 1, 2, and 3. If /z = 0, then Xi = 
constant and z/j = pit + oji, with rji = dFo/Oxi and Ui being constants of integration, and 
a solution is periodic when the values of rji are commensurable; actuall y, if rt = 0 , then 
there are inhnitely many constants Xi that lead to periodic solutions Barrow-Green . 
1997i . 

Then Poincare wanted to know whether the periodic solutions coul d be analytically 


continued when fi remains small. The approach presented below follows [Poincare 


1890 


in which he assumed that the /i 7 ^ 0 solutions at t = 0 are: Xi = ai-\-Sai and yi = Ui-\-Sui, 
with Oj being constants. Now, for t = T, where T is the lowest common multiple of the 
‘In jrji, the solutions are: Xi = ai-\- 6ai -|- Aa^ and yi = Ui-\- ryT -|- 6 ui + Aui, which means 
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that they are periodic if all Aa* = 0 and all Aut = 0. Since only hve equations given 
by Eq. flT^ are independent (as F = const is their integral), Poincare showed that the 
hve equations could be satished if 

where tb is a periodic function with respect to Ui and a; 2 , and if the following additional 
conditions are also satished: the Hessian of d' with respect to uji and UJ 2 must be non-zero 
and the Hessian of Fq with respect to must also be non-zero; the so-called Hessian 
condition is /duiduj 2 ) 7 ^ 0 and similar for Fq. Moreover, if 77 ' = rjiil -|- e:), 

where e is small, then there exists the following periodic solutions Xi = and 

yi = /i, e) with period T' = T/{1 s). 

Using this method, Poincare established that there were periodic orbits for all 
sufficiently small val ues of u. Poincar e’s approach was generalized hrst by the Poincare- 
Bend i xson t heore m Bendixsonl . 


1912 


1901 


, and by the Poincare-Birkhoh theorem Birkhoffl . 
19131. Il915l |. According to these theorems, there are periodic solutions even for 


non-integrable Hamiltonian systems. However, the theorems do not address a problem 
of what happens to solutions of an integrable Hamiltonian system when the governing 
equation s are slightly perturbed? The answer to this question is given by the KAM 
theorem Kolmogorov!. 11954 iMoserl . Il962l. lArnoldl. Il963l | described in Section 3.3. 

An important result concerning t he existence of periodic solutions in the general 


three-body problem was obtained by iHadiidemetrioul [l975a|, who proved that any 


symmetric periodic orbit of the CR3BP could be continued analytically to a periodic 
orbit in the planar general three-body problem Hadjidemetrionl . Il975bl |. Using this 


were constructed 

Bozis and Hac 

iidemetrion. 

1976 

and their stability was investigated 

Hadiidemetriou and Christides. 

I 975 I]. Moreover. 

KatoDodis 

1979 

demonstrated that 


it was possible to generalize Hadjidemetriou’s result from the 3D CR3BP to the 3D 
general three-body problem. 


Periodic orbi 


s in the 3p general three-body problem are typically determined nu- 


mericallv Henon. 

1965, 

1974 

Szebehelv 

1967 

Moore. 

1993, 

Chenciner and Montgomery. 

2000, 

Suvakov and Dmitrasinovid. 

2013 

. To make the computations as efficient as pos- 


sible, the 3D general three -body problem was f ormulated suitably for numerical corn - 
putations [Markellosl . Il98nl | by using an idea of iHadjidemetrion and ChristidesI [1975 . 


Periodic solutions with binary collisions in the general three-body problem were also 
determined analytically for differe nt masses, and 8 d ifferent periodic orbits in a rotating 
frame of reference were found by iDelibaltasI 1983l | . By starting with periodic orbits 


for the 3D CR3BP, families of periodic orbits in the 3D general three-body problem 
were numerically determined, and it was shown that the latter was not isolated but 
instead they for in conti nuous mono-parametric families for given masses of the three 


bodies Markellosl . Il98l| . Let us also point out that there are methods for constructing 


analytically periodic solutions to the general three-body problem by using computer 
experiments 


Valtonen et ah. 1989 
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Figure 3. New periodic solutions are presented on the shape-space sphere with its 
back sides also visible, (a) The figure-8 orbit; (b) Butterfly orbit; (c) Moth orbit; (d) 
Yarn orbit; (e) Yin-yang orbit; (f) The yin-yang orbit in real space. Reproduced after 


Suvakov and Dmitrasinovid 2013j . Copyright 2013 Physical Review Letters. 


More recently, new periodic solutions to the three-body probler n were found 
using combinat ions of analytical and numerical m ethods. Specihcally, iMoord 1993 


discovered, and IChenciner and Montgomery! 2000| independently re-discovered, the so- 
called 8-type periodic solution for the general thre e-body probl e m wit h equal masses 


(see Fig. [3^), and formally proved its existence Montgomery! . 


2001 


A search for 


pe riodic orb i ts in t he vicinity of the 8-type periodic s olution was 
by Suvakov 2013 


Derformed numerically 


The KAM theorem was used by jSimdj [2002l | to establish stability 


and a rotating fra me of reference was adapted to study properties 


of this solution 
of the solution 

soluti ons with some hxed winding numbers and masses were found by IZhang et al. 


Broncke et ahl. I200 (tI |. Moreover, new symmetric non-colli sion periodic 


2004 


Currently, 13 new distinct periodic orbits were found bv ISuvakov and Dmitrasinovic 


2013l |. who considered the planar general three-body problem with equal masses and 


zero angular momentum (see Fig. [3b,c,d,e,f). The authors presented a new classihcation 
of periodic solutions but admitted that their results could not be verihed observationally 
because so far no astronomical systems with the considered properties of the three bodies 
are known. 
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3.3. The KAM theorem and stability of periodic solutions 
Let us write Eq. (fTTll in the following general form 


15 


dq dH dp dH 


(14) 


where q = qi, ..., qn, P = Pi, ■■■, Pn, where t is time. It means that the equations describe 
the autonomous Hamiltonian system with n degrees of freedom. The Hamiltonian is 
H{p, q) = Hoip) + pT[i{p, q) + ..., with H being periodic in q with period 27r, and p is 
a small parameter. 

For the unperturbed motion p = 0, the above equations reduce to 


dq dHo A n 


(16) 


with H = Hi, ..., Qn- Clearly, Eq. flT^ can be integrated where the resulting trajectories 
are conhned to a t orus in th e phase space, actually, t here are invariant tori p = const; 


see 


Barrow-Green 1997 and Diacuand Holmes 1996 for details. For incommensurable 


frequencies the motion is quasi-periodic. Moreover, the system is nondegenerate because 
the Hessian determinant is not zero. 

Now, the KAM theorem tells us that when the system is slightly perturbed most 
of the invariant tori are not destroyed but only slightly shifted in the phase space. This 
has important implications on stability of orbits in the gener al an d restri c ted t hree- 
body problem. The proof of the KAM theorem by Moser 1962 | and Arnold jl963| also 
demonstrated that convergent power series solutions exist for the three-body (as well as 
for the n-body) problem. The KAM theore m seems to be very useful for studying the 
globa l stability in the three-body problem jRobutell. Il993ai iMontgomervl . 1200 ll. [Simo, 


2002 |; however, some of its applications are limited only to small masses of the third 


body and as a result , different meth ods have been developed to deal with the general 
three-body problem jRobutell . Il993b . 


3 . 4 . Collision singularities and regularization 

The Euler, Lagrange and other periodic solutions require suitable initial conditions, so 
an interesting question to ask is: What happens to the three bodies when the initial 
velocities are set equal to zero? Apparently, in this case all three bodies move toward 
their center of mass and undergo a triple collision at that point in hnite time. At 
the point where the triple collision occurs, the solutions (if they are known) become 
abruptly terminated, which means that we have a triple collision singularity at that 
point. Similarly, we may have a double collision singularity when two bodies collide. In 
the following, we describe collision singularities in the general three-body problem and 
present a method called regularization to remove them. 

Based on the important work by IPainleve 


1896. 1897 , it was established that 


collisions between either two or three bodies are the only singularities in the three-body 
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problem, and that such singularities could be removed by setting certain sets of initial 
conditions. The structure of phase space near the singularities changes and solutions 
(if they exist) do not end at the singularities but come close to them and show strange 
behaviours. The work of P ainleve on s i ngula rities resulting f rom c ollisions between 
two bodies was continued by iLevi-Civital 19031 and iBisconcinil 19061 , with the former 
concentrated on the restricted three-body problem and the latter on the general three- 
body problem. They addressed an important problem of initial conditions that lead to 
collisions, and discussed regularization that allows extending possible solutions beyond 
a singularity. Bisconcini deduced two analytical relationships between initial conditions 
and proved that when the relationships were satished a collision did take place in a hnite 
time. 

who 


The problem of triple collisions was investigated by ISnndma.nl |l9n7l . 11909 


formulated and proved two important theorems. First he proved that the triple collisions 
would occur only if all the constants of angular momentum are zero at the same time. 
Second, he demonstrated that if all three bodies collide at one point in space they move 
in the plane of their common center of gravity, and as they approach the collision they 
also approach asymptotically one of the central conhgurations, namely, either Euler’s 
collinear conhguration or Lagrange’s equilateral tri angle. Proofs of Sundma n’s theorems 


and t heir applications are described in detail by ISieeel and Moserl [l99l|. ISundman 


1912 


also studied binary collisions and used the results of those studies to hnd a 
complete solutio n to the general three-body problem (see Section 3.3). Moreover, 


Waldvogell 198nl | considered the variational equations of the three-body problem with 


triple collisions and obtained solutions given in terms of hypergeometric functions, which 
are valid in the vicinity of the collision. 

Regularization means that the motion is extended beyond a singularity, which 
occurs in a solution, through an elastic bounce and without any loss or gain of energy. 
Thus, it is important to know whether such a solution with a singularity can always be 
continued in a meaningful way? ISiegell 194l| showed that in general this cannot be done 
by proving a theorem, which states that an analytical solution that goes through the 
triple collision cannot be foun d for pract i cally a ll masses of the bodies. The work of Siegel 
was signihcantly extended bv iMcGehed 1974| . who obtained new transformations, now 
known as McGehee’s transformations. He also introduced a new coordinate system that 
allowed him to magnify the triple collision singularity and inspect it more closely. What 
he found was a distorted sphere with four horn s extending to inhnity and he named this 


surface a collision manifold; for more details see lDiacu and Holmes 


1996 


There are two 


cu rrently used regularization scheme s, the so-called K-S regularization scheme developed 


bv iKustaanheimo and Stiefel 


Burdet 1967 and 


Heggi^ 


1 


1965 


1973 


and the B-H regularization scheme developed by 


1976|. These two analytical schemes describe close 


two-body encounters and they are typically used to supplement numerical simulations 


Valtonen and Karttnnenl. 


2006 


In Section 1, we briefly discussed the numerical work by ISzebehely and Peters 


1967alJb| . who considered three bodies with their masses proportional to 3, 4 and 5 
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Figure 4. Numerical solutions to the Pythagorean triangle problem show the evolution 
of the system for the first 10 orbits (left) and up to the ejection of the smallest mass 
from the system (right). 


located at the vertices of a Pythagorean triangle with the sides equal to 3, 4 and 5 
length units, and showed that two of these objects form a binary system but the third 
one is expelled from the system (see Fig. 4). Such a pure numerical resul t required 
an an alytical verihcation and a formal mathematical proof that was given by iMcGehee 


1974| for the three-body problem in which the bodies are restricted to move only along 
a hxed line. A more general proof for the planar three-body problem was given by 


WaldvQgell 1976j . The main result of both proofs was that after a close triple encounter 


one of the bodies of the three-body problem would escape from the system. Moreover, 


sufficient conditions for the 

escape of one bodv from t 

le general three-bodv svstem can 

also be soecihed Standish. 

1971. 

Anosova and Orlov. 

1994 

as well as the conditions 

that an ejected bodv would return to the svstem Standish. 

1972 



3.5. Sundman’s complete solutions 


Painleva |1896l . Il897j formulated a conjecture that once collisions are excluded by 


choosing certain initial conditions, then the equations of motion of the general three- 
body problem could be integrable using power series solutions. Painleve neither proved 
his conjecture nor found such solutions. Nevertheless, his results had strong impact on 
his contemporaries who searche d extensively for the solutions but failed to hnd them. 
The only one who succeded was ISnndmanl 19071 1 who found a complete solution to the 
general three-body problem. 

As already stated in Section 1, searches for power series solutions had been 
pursued before Painleve formulated his conjecture. Differen t power se r ies so l utions were 


presented and their convergence rates were investigated jPelaunavl. 118671. Ihindstedt 
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1884. iGvldenl. Il8^|. Some imp ortant points related to the convergence of those 

but he also di d not s uccee d in hnding the solution. 


Poincara. 


1892 


series were made 

A breakthrough occurred when ISnndmanI jlhOTI . 119091 . Il912l | studied double collision 
singularities. Let us now briefly describe Sundman’s complete solutions to the general 
three-body problem by following iBarrow-GreenI 1997 . 

Sundman considered a double collision singularity and showed that it could be 
removed by introducing a regularizing variable u dehned in terms of time t as 


u = 



(16) 


where to corresponds to m = 0 for which the system is regular. With the singularity 
occurring at t = ti, Sundman discovered that the coordinates of the three bodies could 
be expanded in powers of {t — and that he could use an analytical continuation 

to demonstrate that the expansion described correctly motions after the collision. The 
energy remained unchanged, the areal velocity stayed constant, and the solutions were 
valid for t > ti and described motions correctly after each new double collision; note 
that the triple collisions were not allowed in this approach. 

Since Sundman’s regularization transformation was dependent on to, he introduced 
another variable, namely, dt = Tdu, where P = (1 — and 

ro, ri and r 2 are the mutual distances, with the two greater of them being larger than 
I, which gives 0 < P < 1, so there is a one-to-one correspondence between t and u. 

If u* is a real and hnite value of u, then the coordinates of the three bodies, their 
mutual distances and time can be expanded in power series with respect to {u — u*), 
with the radius of convergence G satisfying |a; — ci;*| < G. Finally, Sundman dehned 


r = 


g7ra;/20 _j_ ^ 
g7ra;/2n _j_ ^ ’ 


and used 


2VL 1 + t 

cu = —log-- 

TT 1 — T 


(17) 


(18) 


which allowed him to transform the ca—plane into a circle of unit radius in the r—plane. 
This guarantees that the coordinates of the bodies, their mutual distances and time 
are analytic functions of r everywhere inside the circle. Moreover, the expansions are 
convergent and their different terms can be calculated once I and G are specihed. Thus, 
these are complete solutions to the general three-body problem. 

Glearly, thi s is a remarka.ble result and t his fa. ct had been recognized by Sundman’s 
contemporaries Levi-Givital . Il918l . iBirkhoHi . Il920|. However, there is a major problem 


with Sundman’s solutions; namely, their convergence is extremely slow, actually so slow 
that it requires millions of terms to hnd the motion of one body for insignihcantly short 
durations of time. Since the solutions are also not useful for numerical computations, 
owing to the round-off errors, they have no practical applications. Thus, we have a 
paradox: though we know the complete solutions, nothing is added to the previously 
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accumulated knowledge about the three-body problem! A more detailed discussion 
of Sundman’s complet e sol utions and hi s othe r relevant results can be found in 
Siegel and Moser! 199l| and IB arrow- Greed 1997 . 


3.6. Zero velocity hypersurfaces 


As already mentioned in Section 1 , iHilll 1877L Il878| introduced the zero velocity curves 
(ZVC) for the restricted three-body problem; these curves are boundaries of regions in 
space where the bodies are allowed to move. Hill’s results can be generalized to the 3D 
general three-body problem; however, in this case, regions of space where motions are 
disallowed are not curves but rather 4D surfaces, or hyper surfaces. The so- called zero 
veloc ity hypersurfaces (ZVH) can be determined using Sundman’s inequality jSergyselsI . 
1986| and the resulting restrictions on the motions must obviously be consis t ent w ith 


the existing solutions to the problem Marchall. 


1990 


Valtonen and Karttuned . l2006| . If 
the three bodies have masses Mi > M 2 > M 3 , and if the terms of the order {M^/M 2 Y 
can be neglected, then the ZVH of the general three-body p roblem become the ZVC of 
the restricted three-body problem Milan! and Nobili . Il983 |. 


4. Numerical solutions to the general three-body problem 


Since the analytical solutions to the three-body problem obtained by ISnndman 


1907 


converge extremely slowly, the modern practice of solving this problem 
involves numerical computations. In general, there are different numerical methods 
to solve the three-body problem, and determine periodic orbits, resonances and 
chaotic indicators. In the following, we present an overview of these methods, 
including numerical simulation packages such as MERCURY, SWIFT, and HNBody 
developed for astronomical applications Chambersl . Il999l. iLevison and Dnnca.nl . 11994 


Ranch and Hamilton! . 12002 


4 . 1 . Numerical methods 

One of the commonly used methods in determining numerical solutions is the application 
of symplectic integration. This method relies on the mathematical parametrization of 
the equations of motion into matrices and assumes a low occurrence of events where 
the orbital velocity could quickly become large. In celestial mechanics, this is an 
implicit assumption on the likely values of eccentricity for the system being considered. 
With low eccentricity, larger time steps can be implemented, which can signihcantly 
quicken the calculations. This is possible due to the symplectic nature of Hamiltonian 
dynamics where matrices are symmetric and easily invertible. Using this mathematical 
property, one can study nonl inear Hamiltonian dynamic s with computationally fast 
mapping integration methods. IWisdom and HolmanI 1992l | developed an algorithm that 
incorporates this idea and applied it to study the n-body problem. 
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In general, one can envision a basic mapping method such as Euler’s or leapfrog’s 
integration, implemented with symplectic properties to attain higher speed and accuracy 
over more compu tationally expensive methods such as the Runge-Kutta method 


Press et al.l. Il992| . We now briefly outline the leapfrog method, and then demonstrate 


how it can be transformed into a symplectic method. 

The leapfrog integration method applies to problems that can be expressed as a 
second order differential equation of the following form: 


X = F{x). (19) 

Many celestial mechanics problems are of this form and could thus be easily solved 
using the leapfrog integration method, which is stable for oscillatory motion as long 
as the time step. At, is constant and less than 2/a;; this parameter is related to the 
mean motion u = ^JG{M + m)/a^, where a is a semi-major axis. For example, the 
numerical integration of Solar System dynamics has a maximum time step on the order 
of 0.0767 years. However, in practice, one will use a smaller time step to better sample 
Mercury’s orbit. The choice of At will be at least an order of magnitude larger than 
other numerical methods and hence computationally quicker. The leapfrog method of 
integration utilizes half step evaluations of the following form: 


Xi = Xi-i + Vi-i/2At, 
tti = F{xi), 

Vi+l/2 = Vi-i/2 + ttiAt, (20) 


where Oj are coefficients determined by F[xi), and Vi is the velocity. 

The main advantage of this method is its symplectic properties, which inherently 
conserve the orbital energy of the system. This is highly valued in the held of celestial 
mechanics. Also the leapfrog method is time-reversible, which means that one can 
integrate forward n steps to a future state and reverse the direction of integration to 
return to the initial state. 

Other commonly used computational techniques can be categorized broadly as 
Runge-Kutta like methods. Different approaches to the general algorithm of Runge- 
Kutta have produced methods of determining the coefficients for the Runge-Kutta 
method to high order but with a decrease in computational speed. Such approaches 


inclu de Bulirsch-Stoer, Dormand-Prince, and Stormer-Cowell techniques Press et ah 


1992 


been 


where variable time steps along with sophisticated predictor-correctors have 
utilized to make them more adaptive and accurate when compared with 


symplectic methods. We show the basic approach of the Runge-Kutta-Fehlberg (RKF) 
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implementation following iPress et al.l 1992| . The basic generalization is 


Vn+l Vn T ^ ^ bi ki, 
i=l 

where 

ki = h f(tn,yn) 

k2 = h f{tn + Cah, yn + Oai/ci) 


ks — h f (tn + Cgh^yn + Clslfcl + Cis2k2 + • ■ • + ttg^s-lks-l), (21) 


where yn is a dependent variable, is an independent variable, s is a number of 
stages, and the coefficients bi, Ci and Oij, with i = l,2 ,...,s and 1 < j < i < s, are 
determined from the so-called Butcher tableau, which give the relationships between 
these coefficients. The most common RKF implementation is RKF45, which describes 
a method with s = 6. The adaptive portion of the method compares an order m method 
{yn+i) with an order m — 1 approach where m = 5 for RKF45. This comparison 

is used to determine the relative error between the orders, which will be a criterion 
to modify the step-size h. The error estimate e„+i with the control scalar ^ is of the 
following form 


Cn+l Vn+l Vn+l 


h eh Y 

V2|en+iiy 


i=l 


( 22 ) 


In order to determine the value of z, one specihes the error control tolerance e. 
Then, the product of 2 ; with the step-size h is used as the new optimal step-size h' 
for the next iteration. When considering simulations where close encounters are highly 
probable, an adaptive step-size method is preferred. Otherwise large errors in the total 
energy can accrue and impact the reliability of the results. 

A more classical method of integration uses the Taylor method, which relies on the 
knowledge of terms within a Taylor series. In general, the optimal number of terms in 
the Taylor series may not be known, thus optimization packages such as TAYLOR uses a 
text input of the ODEs to be solved and prov ides the user with an optimal time-stepper 


that will minimize error Jorba and Zod . 120051 . This method can be comparable in terms 


of speed and accuracy with the other methods presented; however implicit assumptions 
are made that would need to be evaluated before a specihc application. 


4-2. Numerical search for periodic orbits and resonances 

Periodic orbits have been determined analytically for special cases with the Euler and 
Lagrange solutions being the most notable. However, other analytical solutions have also 
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been found in other periodic families by using numerical methods (see ISuvakovi 2013 
for de tails) . Studies of periodic orbits are related to the topic of resonance. iMardling 
j2nn8j | and Murray and Dermott 1999 | give detailed descriptions of the perturbation 
theory involved in determining the locations of resonance based on approximations to 
a harmonic oscillator. Using the perturbation theory, appropriate initial guesses can be 
made as to the locations of periodic orbits and resonances within a parameter space 
of semi-major axis and eccentricity. Such guesses are made using the ratio of Kepler’s 
harmonic law between two orbiting bodies. 

The nominal resonance location, Ra, can be dehned in terms of two integers {k, /), 
where I is the order the resonance to form a ratio of k/{k 1). For internal resonances, 
the third body orbits between the primary and secondary, and the relation is given by 


Ra = 


k 


k + l 


2/3 


rUr, 


rUj, -|- m. 


(23) 


where a, rrip, and denote the semi-major axis of the secondary, mass of the primary, 
and mass of the secondary, respectively. There is a similar dehnition for the external 
resonances. 

Knowing this relation, one can begin to search for resonances that lead to periodic 
orbits using the general form of the resonant argument as 


(j) — jiX' j2\ + 5 


(24) 


with conditions on the integer coefficients using the d’Almebert rules Hamiltoi] . 


1994| . which require a zero sum of the coefficients and an even sum for the and 
je coefficients due to symmetry conditions. Also, the A', w', 12' represent the mean 
longitude, longitude of pericenter, and ascending node of the perturber. Accounting for 
all these conditions, a numerical search for resonances can be performed in a tractable 
way using the period ratio P'/P as an initial guess. 

Eq. [21] can be characterized in terms of mean motion and secular components. 
When considering mean motion resonances, the ji — coefficients are included such 
that they obey the zero sum condition, and the remaining terms are assumed to be 
zero. The search for secular resonances usually consider only the condition where 
js = —ji = 1 and all other coefficients are zero. The differences between these types 
of resonances manifests in the changes of the Keplerian o rbital elements, which dehn e 
an orbit. Further details on these changes can be found in [Murray and Dermott! 1999l |. 
and their applicability of study within the Solar System. 

Another method of searching for periodic orbits includes the characterization of 
chaos in the three-body problem. Chaotic regions can alter the periodicity of the orbit, as 


Gozdziewski et ah 

2013 

and 

Misraszewski et ah 

12912 


this principle with a chaos indicator (for details, see Section 4.5). 



































The three-body problem 

4-3. Maximum Lyapunov exponent 


23 


Time evolution of dynami cal systems can b e characterized through the use of the method 


of Lyapunov exponents jLvapnnovl . Il907|. Firs t applications o f these exp onents to 
the general three-body problem were made by iBenettin et ahl 19761. Il980|: specihc 


applic ations to the restricted three-body problem were hrst made by iJeffervs and Yi 


1983 


The method of Lyapunov exponents utilizes a description of two nearby 
trajectories in phase space for a third mass within the phase space. The mathematical 
description follows a power law (usually linear) for stable orbits, whereas an exponential 
law is used for unstable orbits, with the Lyapunov exponent being a multiplicative power 
of such laws (i.e., 

Most generally, the rate of divergence resulting from these laws is divided into 
3 regimes. A positive Lyapnuov exponent A denotes that a trajectory, corresponding 
to one degree of freedom, is diverging from the neighboring orbit, with the rate of 
divergence being characterized by a time series A(f). A negative Lyapunov exponent 
implies that a trajectory is converging onto a stable manifold and can also be 
further described by a time series. The third case considers a Lyapunov exponent 
that is exactly equal to zero, which means that the two trajectories are parallel to 
each other and by induction will remain so unless another force ar ises to disrupt 
the system. The mathemati cal dehnition of the Lyapunov exponent Hilbornl . Il994l . 


Musielak and Musielakl . 


2009 


IS 


A = lim 


M. 

t^oa t — to \5{to) 


In 


(25) 


Using the above basic description, one can calculate the spectrum or set of Lyapunov 
exponents for a given system with n degrees of freedom in a 2n dimensional phase space. 
Within the CR3BP, the system can be rotated with respect to the motion of the two 
large masses, so n = 3 degrees of freedom exist within a 6 dimensional phase space. 
Since in this example, the system is Hamiltonian, the trace of the Jacobian J must be 
zero, which presents a symmetry among the Lyapunov exponents, with 3 positive and 
3 negative exponents existing. There are other examples for which such symmetry does 
not exist, and only 2 regimes, chaotic or dissipative, can be identihed. The sum of the 
exponents (Tr J) can then be positive or negative, and this sum determines whether a 
system is chaotic or dissipative, respectively. 

The Jacobian for the Ah mass is given by 


J^ 


( 0 
0 
0 

dji 

or], 

\di 


0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

dji 

Ml 

0 

0 

dvi 

dQ 

di}i 

dvi 

0 

0 

diji 

dii 

dQ 

dCi 

0 

0 

drji 

dCi 


0\ 

0 

1 

0 

0 




(26) 
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where the coordinates {^,ri,() denote an inertial reference frame, 
are given by Newton’s Universal Law of Gravitation, i.e., a* = 


The accelerations 


-E 


GMi 




and 


4■4- Feist Lyapunov Indicator 


Since the calcnlation of the Jacobian involves extra computational time, other methods 
have been developed to determine more efficiently the onset of chaos. We discuss two 
such methods: the Fast Lyapunov Indicator (FLI) and the Mean Exponential Growth 
Factor of Nearby Orbits (MEGNO); the latter is discussed in the next section. Using 
these methods chaos can be detected within 100,000 year timescales; note that previous 
methods required much longer timescales. 

The FLI is a method, which identihes two shortcomings of the previous 
approach of using the Lyapunov exponents, namely, normalization and dimensionality 


Froeschle et ahl. Il997l |. The problem of dimensionality is related to the computational 


time used to follow non-chaotic vectors within the phase space. Thus, the supremum of 
the tangent vectors is used to determine the chaos indicator, with other vectors being 
neglected. This is shown by the following standard formula: 


FLI(t) = supj , (27) 

where j = 1... m of a m-dimensional basis. 

The slope of the time series of FLI are used as indicators of chaos. If the slope 
is steep and positive, then chaos is inferred; this is analogous to a positive Lyapunov 
exponent. The slope may also remain flat, neither increasing or decreasing, which is 
the stable case. Lastly, there is the possibility for a negative slope, but this is not 
generally realized in celestial mechanics as the systems under study are Hamiltonian 
and dissipation would imply a loss of energy. 


4-5. Mean exponential growth factor of nearby orbits 


An alternate technique for identiheation of chaos is the use of the MEGNO; this method 
is based on the maximum Lyapunov exponent but in a sli ghtly different manner. The 
method was developed by [Gincotta and Simdl . Il999l. l2000| as a general tool in probing 
chaos of Hamiltonian systems. Since t hat time, it has been used to probe the dynamics of 


Lauehlin et ah. 

2002. 

Pilat-Lohineer et ah. 

2 OO 3 I . For exan 

iples 

using the MI 

IGNO 

met ho ( 

; and other chaos indicator methods, we recommend 

Pavlov and Macieiewski 

2003 . 

Kotoulas and Vovatzis 

2004 


Barnes and Greenberg 

2006 

Funk et ahl 

2009 . 


Dvorak et al.l 


2010 


Hinseetal 


2010 


Mestre et ah 


2011 


, and ISatval et al 


2013 


The MEGNO method employs the integral form of determining the maximum 
Lyapunov exponent. In this sense, it can be interpreted as a measurement of the time- 
average mean value of the hnite time Lyapunov exponent with the following dehnition 
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Gozdziewski et al.l. 120011. IHinse et al.l. 12010 


Y{t) = 


along with its time-averaged mean value 


h(s) 


s ds, 


(28) 


(r)(«) = l f‘Y(s)ds. 


(29) 


The advantage of this representation is that (Y) (t) co nverges faster to its limit value 


even with the cost of two additional equations of motion Gozdziewski et al.l . 1200 1| . This 


is due to the fact that the time-weighting factor amplihes the presence of stochastic 
behavior, which allows the early detection (i.e., shorter integration timescales) of chaos. 
The additional differential equations are simply the differential forms of Eqs. fl28|) and 
given respectively as 


6 X 

X = -t and w = 2—, 
h f ’ 


(30) 


where Y{t) = 2x(t)/t and (Y) (t) = w(t)/t, with all the variables having their usual 
meaning. 


5. Applications of the general three-body problem 

The general three-body problem has been considered in many astronomical and 
spaceflight settings, including galactic dynamics, stellar formation, and also the 
determination of trajectories for spacecraft missions, such as manned, satellite, and 
robotic landers. 


5.1. Astronomical settings 


Exploring the general three body problem within an astronomical setting can encompass 
a broad scale. Here we discuss only the most basic cases. In order for a three-body 
problem to be considered general, all three masses must interact with each other with 
no limitations on the eccentricity of their orbits. One example is the case of triple 
stellar systems that are not hierarchical. The basic dynamics of these systems can 
only be investigated numerically because a variable timestep method of integration (see 
Section 4) is required due to Kepler’s Second Law. The masses in such settings are likely 
to undergo several close approaches and great care must be taken during the point of 
close encounter between any of the masses. 

As previously discussed (see Section 3.4), the Pythagorean problem exemplihes 
the difficulty in obtaining meanin gful results in the face of lar ge sources of numerical 
error. In terms of stellar dynamics, iReipnrth and Mikkolal 2ni2| investigated this setup. 
Examples include a hnite sized cloud core around the three masses, which presents an 
even more intricate dynamical problem as the masses of the three bodies increase with 
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time depending on the interaction with the cloud core. These types of investigations 
provide an insight into the formation of wide binary stars (> 100 AU separations) due 
to dynamical instabilities in the triple systems. 

Conversely the general three body problem has been used to explain the 
formation of close binaries (< 100 AU separations). Results of numerical evaluation 
of orbits along with secular perturbat ions have explained the occurrence of close 


binaries Fabrvckv and Tremainel . l2007l |. who provided a specih c mechanism for the 


form ation of such systems through the Lidov-Kozai cycles Lidovl . Il962l . iKozail . 


1962 


The 


; within these cycles large variations of eccentricity and inclination occur, 
idea that stellar binary system may host planets has been explore 


m 


great detail iGonczi and Froeschle. 1981, Rabl and Dvorak 

00 

00 

rH 

Dvorak et ah. 

1989, 

Holman and Wiegert 

1999. 

Pilat-Lohinger and Dvorak. 

IOO 2 I. 

Pilat-Lohinger et ah. 

20031. Dvorak et ah. 

2004, 

Vlnsielak et ah. 2005. Pilat-hohinger and Dvorak. 

2008, 

Haghighinour et al.. / 

1010, Guntz.|2014|. with the additional observational verihcation bv 


the NASA Kepler mission (see Section 1). The consideration of binaries hosting planets 
will be addressed in more detail in Sections 8 within different regimes of approximation. 


5.2. Exoplanets 

Recently, using the Kepler data, over 700 exoplanet candidates have obtained the status 
of exoplanets based on a stat istical method to identify the likelihoo d of false positive 


candidates among the sample (Lissaner et al.l. 12014 iRowe et al.l . l2014| . This portends to 


the general three body problem as each considered case includes the possible interaction 
of three masses. This is due to the possible geometric conhgurations in the photometric 
method. Some examples include the presence of third light contamination, background 
eclipsing binaries, and stellar triple systems. Prior to this major identihcation of 
exoplanets, a method called transit timing variations (TTVs) has been used to identify 
systems with multiple planets (Kepler-11, Kepler-36, Kepler-9); the method determines 
how three- body interactions would affect th e timing of transits r elativ e to a linear 


prediction [Lissaner et al.l. l201lL ICarter et al.l. l2012l. iHolman et al.l. |2010|. We presen t 


the basics of TTVs and defer the reader to a more in depth description by IWinnI 2011 


The method of TTVs relies on the observable deviation from Keplerian dynamics. 
The deviation reveals the existence of additional bodies, which in principle is a similar 
method that was used in the discovery of Neptune throug h variations in the orbit of 


Uranus AdamsI . ll846L lAirvi . ll846L IChallisl. Il846l. iGalM Il846| . Using the times of transit 


and the planetary period of a transiting planet, a htting is made to a linear function to 
predict the subsequent eclipses (f„ = to + nP) of that planet. When the planet eclipses 
a host star early or late relative to the predicted time of the htted funct i on, a TTV 
is observed. The power of this method was illustrated by iNesvorny et al.l 2ni2| with 
the detection of a non-transiting planet (KOI-872c) inferred by observations of the seen 
transiting planet (KOI-872b). The vetting of this planet underwent intense scrutiny 
from the perspective of the three-body problem as only particular solutions allowed for 
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the determination of the planet’s mass as well as keeping the perturbing planet stable. 

Many other discoveries have been made detecting multiple planet systems, including 
KOI-142 which i s anot her discovery of a planet pair based on their mutual interactions 


Nesvornv et al.l. |2013| . Much previous work was related to predicting the degree to 


which discoveries could be made based on the possible pertu rbations, and involved 
determining the sen si tivity required to perfor m observations Miralda-Escuda . 


Holman and Murrav. 

2005. 

Agol et ah. 

2005 


2002 


TTVs along with transit duration 
variations (TDVs) ha ve even b e en su ggested as possible methods to detect the presence 


of extrasolar moons Kippingl. |2009| . although none so far have yet to be observed. 


Additional discussion on the requirements for the detection of extrasolar moons is 
presented in 1 18.51 

To account for exoplanets and exomoons in extra-solar planetary systems with 
the number of objects higher than three, additional new tools, such as TTVfast and 


yhotodynam, have been deve l oped with full n-body capabilities built-in [Deck et ah 


2Q14, 

Carter et al.. 

2Q11, 

Pal, 

2012 


5.3. Spacecraft trajectories 

Spaceflight applications involve the determination of spacecraft trajectories for satellite 
and lander space missions. NASA’s Apollo 8 - id missions to the Moon required accurate 
numerical solutions to the three-body problem. NASA also launched several discovery 
missions, with the most notable being the Hubble Space Telescope, CHANDRA X-Ray 
Observatory, SPITZER Space Telescope, Kepler Space Telescope, and the upcoming 
James Webb Space Telescope. Each of these missions involve different solutions to three- 
body problem and all but Kepler orbit the Earth. One of the most interesting satellites 
will be the Transiting Exoplanet Survey Satellite (TES^) as it will be in a high-Earth 
orbit and in a 2:1 resonance with the Moon. For these applications, it is required 
to consider the motions of the terrestrial objects relative to each other, including the 
motions of the Moon and Mars. The orbits of both of these objects can be affected by 
the perturbations from other planets. 

More basic applications of the general three body problem involve artihcial satellites 
commonly used for communication and/or military purposes. These satellites are 
generally considered within the realm of a restricted three-body problem (see Sections 
6, 7 and 8). A restricted three body problem simplihes the solution because one of the 
masses is much less than the others and has negligible gravitational influence. Thus, 
the much larger masses have approximately Keplerian orbits. This is the case for the 
Earth-Moon-satellite system, and it is discussed in detail in Section 8. 


II http://tess.mit.edu/ 
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Figure 5. Non-r otating and rotatini 


Reproduced after lEberle et al 


200 


coordinate systems used to describe the CR3BP. 
Copyright 2008 Astronomy & Astrophysics. 


6. The circular restricted three-body problem 


6.1. Governing equations 

Based on the criteria stated in Section 1, the circular restricted three-body problem 
(CR3BP) requires that the two objects have their masses signihcantly larger than the 
third one (Mi >> M 3 and M 2 >> M 3 ), and that the motions of Mi and M 2 are 
limited to circular orbits around their center of mass. In the literature devoted to the 


CR3BP, Ml and M 2 are typically referred to as the primaries jSzebehelvI . Il967l |. but in 
stellar dynamics Mi and M 2 are called the primary and secondary, respectively, with 
the assumption that Mi > M 2 . In this paper, we use the former nomenclature when 
a clear distinction between Mi and M 2 is not necessary; however, we use the second 
nomenclature when a distinction between Mi and M 2 becomes important in describing 
a specihc astronomical application. 

Now, if the third mass moves in 3D, this case is called the 3D CR3BP; if it moves 
in the same plane as the primaries, we call this case the planar CR3BP. Since the 
gravitational influence of the third body on the primaries is negligible, the orbits of 
the primaries are described by the two-body problem, whose solutions are well-known. 
Once the solutions for the primaries are known, they can be used to determine the 
motion of the third body resulting from the gravitational held of the primaries. From a 
mathematical point of view, the problem becomes simpler than the general three-body 
problem. Nevertheless, for some sets of initial conditions the resulting motions of the 
third body remain unpredictable. 

Let us choose the origin of a non-rotating coordinate system (see Fig. ED to be at 
the center of mass of the primaries, so we can write MiRi = M 2 R 2 , where Mi and M 2 
are the masses of the primaries with Ri and R 2 denoting their position vectors with 
respect to the origin of the coordinate sys tem. To describ e the problem mathematically, 
we follow Szebehelv 1967], Danby 1988 | and Rov 20^. 

The components of the position vectors are: Ri = (Xi,y'i,Zi), R 2 = {X 2 JY 2 , Z 2 ) 
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and R 3 = (XsjYs, Z 3 ), and they are nsed to express r^i = Ri — = ri and 

^32 = -R 2 ~ -R 3 = h *2 in terms of their components. The eqnations of motion for the 
CR3BP can be written as 

^ . M: ^ 

(31) 


<eR, 


dR 




‘(R,-Rz ), 


i=i 


where r, = [{X, - + (Yj - F 3 )' + {Zj - Z,] 


2]l/2_ 


We dehne M = Mi + M 2 , /i = M 2 /M and a = 1 — n, which give Mi = aM and 
M 2 = fiM. With the gravitational force being eqnal to the centripetal force, we have 


= Ry 


(32) 


where D = Ri + R 2 , u is the angnlar freqnency or mean motion, and j = 1 and 2. Since 
M 3 _j = MRj/D, Kepler’s third law = GM/D^ is antomatically obtained. Moreover, 
we also have Ri = nD and R 2 = aD. 

Now, once the large two masses are in circnlar orbits, we write: Wi(f) = joD cos cat, 
Yi(t) = piD smojt, X 2 (t) = —aD cosojt and Y 2 (t) = —aD sincnf, with Zi = Z 2 = 0, if 
the orbits are in the same plane. 

Introdncing X 3 = Dx, Y 3 = Dy and Z^ = Dz, and allowing the third body to move 
in 3D, we nse Eq. flMj) to derive the following set of eqnations of motion 


.. a , . u , . 

X = — Dx — ficosT) - Dx + a cost) , 


a , . 3 h / -3 

y = — Dy - /ismr) Dy + asmr) 


and 


z = 


a 

Q T Q 


(33) 

(34) 

(35) 


where r = ojt, and x, y and z represent the second-derivative of x, y and z with 
respect to r, respectively. This set of eqnations describes the CR3BP in the non-rotating 
coordinate system. 

We write the above set of eqnations in the rotating (synodic) coordinate system 
by nsing the following relationships between the coordinates x, y and 2 ; in the inertial 
system, and the coordinates x*, y* and z* in the synodic system: x = x* cos t — y* sin r, 
y = X* sinr-|-?/* cost and = z*. Therefore, the set of eqna t ions d escribing the CR3BP 
in the synodic coordinate system can be written as Eberlel . 


2010 


X 


2y* = X* 


a 


{x* - fi) - + a) , 


y* + 2x* = (1-^- 


a 




y 


(36) 

(37) 


and 


z* = 


a 

Q H“ Q 


( 38 ) 


where ri = D[{x* — /i)^ -|- {y*Y + (z*)^]^/^, and r2 = D[{x* -f a)^ -f {y*Y -|- (z*)^]^/^. 
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Lagrange 1772| found interesting solutions to the CR3BP that describe equilibrium 
positions of the third body when all net forces acting on it are zero. In this case, Eqs 
, iiD, and fl38|) reduce to 


X* — — fi) — a) = 0 , 


(39) 


a 


r 1 - - 




= 0 


and 


a 


+ 




= 0 


(40) 


(41) 


1 ' 2 . 

Taking 7 ^ 0 gives y* = 0; however, z* = 0 gives y* 7 ^ 0, which means that the 
equilibrium solutions are conhned to a plane. We take this plane to be the a;*j/*-plane. 
Thus, with y* 7 ^ 0, we obtain ri = r 2 = 1, which shows that the third body can be at 
either of these two points loc ated at the same distance from the primaries as the two 


primaries are from each other Rovl . |2005| . Clearly, an equilateral triangle is formed and 


the two equilibrium points are called the Lagrange triangular L 4 and L 5 points (see Fig. 
6 ). Both points are stable for the mass ratios 0 < /i < po, where /iq = (1 — •\/^/9)/2 is 
called Routh’s critical mass ratio Bard ini . |2002| . Now, with ?/* = z* = 0, it is easy to 
show that Eq. fl3^ has three equilibrium solutions located on the line passing through 
the primaries. These three solutions are called the Lagrange collinear Li , a nd L 3 
points (see Fig. 6 ), and they are unstable for any value of /i jRagos et al.l. |2000| . The 
Lagrange points are also called the libration points. 


6.3. Jacobi’s integral and constant 

Let us introduce a potential function 0* dehned as 


= i [(V)^ + (!,-)^] + f + f 

/ ri r2 


(42) 


and write Eqs. (136|) . fl37)) and 


in terms of 
X* - 2y* = 
ij* + 2x* = 


and 


2 : = 


d(j) 


dz 


d(j)* 
dx* ’ 

(43) 

d(j)* 

dy* 

(44) 

>* 

* 

(45) 


We multiply Eqs. (03]), fl44D and fl43]l by x*, y* and i*, respectively, sum them up, 
and integrate them to obtain the following Jacobi integral 


(n*)2 = 20* - C* , 


( 46 ) 
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where (f*)^ = {x*Y + {y*Y + and C* is the Jacobi constant in the rotating 

coordinate system given by 


C* — /i + 2/ipo H- - + — - -1" 2-^/po(l ~ h) 


(47) 


Po 1 + Po 

Here po = Rq/D (see Fig. 5) represents the normalized initial position of the third mass. 

Having obtained the Jacobi integral in the rotating coordinate system, we also write 
its explicit form in the non-rotating coordinate system 


— 2{xy — yx) = 2 { — + — \ — C , 

Ji r2 ' 


(48) 


where v'^ = x'^ -\- y"^ -\- and C is the Jacobi constant in the non-rotating coordinate 
system. 

6.4- Zero velocity curves 

In the case when the particle’s velocity v* is zero, the Jacobi integral given by Eq. (1461) 
reduces to 2(f)* = C*, which dehnes the zero velocity curves (ZVCs) for the third body. 
Once C* is determined by the initial conditions, the ZVC encloses a region in space 
where 20* > C* and where the third body must be found. The third body cannot be 
outside of this region because its velocity would be imaginary since 20* < C*. The 
ZVCs do not give any information about orbits of the third body in the region of space 
where it is allowed to move. 

The shape of the ZVCs is uniquely determined by the values of C*. Changes in the 
shape of the ZVCs in the x*y* plane with changing C* expressed in terms of the initial 
position of the third mass po (see Eq. 07]) are plotted in Fig. [6l the region inaccessible 
to the third mass is colored. The presented results were obtained for p = 0.3 and po 
ranging from 0 . 2 , which corresponds to a large value of C*, to 0 . 6 , which corresponds 
to a small value of C*. Note that the regions inaccessible to the third mass shrink and 
eventually vanish at the Lagrange points L 4 and L 5 , as shown in Figs. [ 6 lJ,e,f. 

Transit orbits in the CR3BP are the trajectories that pass through the neck region 
of the ZVCs, which occurs in the vicinity of the Lagrange points Li, L 2 , and L 3 where 
the ZVCs open at these libration points at certain critical values of C* (see Fig. jej). 
Moving along these transit orbits, the third body can transit between the primaries or 
even escape from the system. Extensive numerical studies of transit orbits in the CR3BP 


were recently performed by iRen and ShanI 201^. They established the necessary and 


sufficient conditions for such transition by using manifolds of the vertical and horizontal 
Lyapunov orbits along with the transit cones, and applied the results to the Sun-Earth 
CR3BP. 


6.5. Periodic and quasi-periodic orbits 


Poincare 1892l | attached great importance to periodic orbits in the CR3BP, and he used 
his method (see Section 3.2) to identify and classify periodic motions of the third body 
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Figure 6. Zero velocity curves in the synodic coordinate system, computed for a mass 
ratio y = 0.3 and various initial positions po = Ro/D of the third mass. The dashed 
lines represent the boundary of zero velocity while the coloured region indicates the 
forbidden zones, (a) po = 0.2; (b) po = 0.28; (c) po — 0.36; (d) po = 0.44; (e) po = 
0.52; (f) Po = 0.6; 


with respect to the synodic coordinate system. He searched for periodic orbits by taking 
/i = 0 and then used analytical continuation to hnd periodic orbits for fi > 0. In his 
classihcation of periodic orbits, premiere sorte contains the orbits generated from the 
known two-body circular orbits; deuxieme sorte contains the orbits generated from the 
known two-body elliptical orbits; and troisieme sorte contains the orbits generated from 
the known two-body circular orbits, but with a non-zero inclination of the third body 
with respect to the plane in which the orbits of the primaries ar e conhn ed. _ 


Moulton et ah 

1929 

and 

Stromgren and Pedersen 

1922 

Then, 

Stromgren 

|l935 


performed extensive studies of the periodic orbits for /i = 0.5 and demonstrated that 
they belong to the premiere sorte. He concluded that the termination of this family of 
orbits must be an asymptotic periodic orb it spiraling into the Lagrange poin ts L 4 and 
L 5 ; a proof of this conjecture was given bv iHenrardI 1973| and iBuffoni 


1999 


Various families of periodic orbits in the commensurable 3D CR3BP, where two 
bodies are close to commensura bility in t heir m ean motions around the central body. 


were identihed and classihed by ISinclairl jl970j . Moreover, a special class of periodic 


orbits in the CR3BP resulting from an analytical continuation of Keplerian rectilinear 
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periodic motions was discovered by iKurcheeval |1973j , who showed that the period of 
such orbits is an analytical func tion of u\ the conditions when these orbits are stable 

Comprehensi ve studies of the periodic 


995 


were established analytically by lAhmad 
orbits in the CR3BP were performed by 
also investigated intersections between families of generating orbits Henonl . ll997 . 


Henon 1965. 19741 who classified them and 


29911 


1935 


There are periodic orbits around the Lagrange points as demonstrated bv IStromgren 
Periodic orbits in the vicinity of these points are typically classihed as halo, 


vertical, and horizontal Lyapunov orbits Henrard and Navarrol . l2994j . They arise either 
as a continuation of inhnitesimal oscillations or as bifurcations of the planar periodic 
orbits around the points. Some unstable periodic orbits around Li, L 2 and L 3 were 


identihed and classihed by 

Zagouras an 

Kazantzis 

1979 

Since L 4 and are stable Sz^b^h^ly. 

1967, 

0 

2995 I 


there are stable zones around 


, 1991 

. The existence 

1. Perdios et ah 

1991 


identihed the long and short periodic orb its near the two libration points for the mass 


parameter fi ranging from 9.93 to 9.5, and iMarkellos 


1991 


found families of remarkable 


termination orbits. According to Routh’s criterion, the existence of periodic orbits near 
the triangular libration poi nts is limited to the mas s ratios 9 < u < Po, where po is 


Routh’s critical mass ratio Bardinl. |2992| . Moreover. iBrouckd 1999| adopted canonical 


units to dehne the unit circle and formulated a symmetry theorem that was used to 
discover new long and short-period orbits around L 4 and L 5 . 

All orbits discovered in the above studies were symmetric periodic orbits, which 
are much more easy to generate than non-symmetric periodic orbits. Systematic 
explorations of non-symmetr ic pe riodic orbits in the 3D CR3BP were done by 

who discovered new orbits 


ZagQuras and MarkellosI jl977l | and IZagouras et ah 


1996 


and classihed them. Searches for non-symmetric periodic orbits near L 4 (because of the 


symm etry conditions the orbits are the same at L 5 ) were performed by iHenrardI 1979 


2992| . who discovered an inhnite number of families of non-symme t ric pe riodic orbits 


in the vicinity of the points. In recent work of lHenrard and Navarrol 2994| showed that 
some of these orbits are emanating from homoclinic orbits (see Section 6.7). 

The existence of quasi-periodic orbits around the libration points was conhrmed by 


1999. 

Gomez et ah. 

1999 

. An overview of 

iHenrard and Navarro 

2994 

, who showed 


that different families of such orbits exist and that these families are associated with 
the known halo, vertical and horizontal Lyapunov orbits. They also presented a fast 
numerical method based on multi ple Poincare sections (see S ection 6.7) and used it 
to hnd quasi-halo periodic orbits. 


Henrard and Navarro 29941 demonstrated that the 


method gives full convergence for a given family of quasi-periodic orbits, and that it is 
robust near chaotic regions. 
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In the Hill three-body problem, the mass of one primary dominates over the other 
masses (Mi >> M 2 ), and the mass of the third body is negligible when compared to 
both primaries (Mi >> M 3 and M 2 >> M 3 ). Because of these mass relationships, 
M 2 moves around Mi, and if its orbit is circular, the problem is called the circular 
Hill problem. The third body can move either in 3D space or in the same plane as 
the circular orbit, and its or bit is determined by the equations of motion governing 
the problem Szebehelvi . Il967l | . The equations of motion in the non-rotating coordinate 
system are obtained by taking p — )■ 0 in Eqs (l3^ through ([35]), and in the synodic 
coordinate system by applying the same limit to Eqs fl36l) through fl38l) . 

An extensive discussion of the circul ar Hill problem a nd its specihc applications 
to the Sun-Earth-Moon system is given by ISzebehelyl 19671], who also presents families 
of periodic orbits, the libration points and the ZVCs. Applicatio ns of the circular Hill 
proble m to the Sun-Jupiter-asteroid system was considered by iBolotin and MacKay 
2nnnj . who proved the existence of an inhnite number of periodic and chaotic (almost 
collision) orb its in the synodic coordinate sy stem. Moreover, so me periodic orbits were 


identi hed by IZagonras and MarkellosI 1985| , and iHenonI 20031 1 extended his previous 


work HenonI , Il974| on the circular Hill problem and discovered new families of double 
and triple-periodic symmetric orbits. 


6.7. Poincare’s qualitative methods and chaos 


Poincare 1892l | considered an unstable periodic solution to the CR3BP (or in some 
specihc cases to the circular Hill problem) and represented it by a closed curve. He 
also accounted for a family of asymptotic solutions located on two asymptotic surfaces. 
The closed c urve and the two asym ptotic surfaces were his geometrical representations of 
the CR3BP Barrow- Greenl. Il997l | . Poincare demonstrated that the equations describing 
these surfaces could be written as 


— = f(yi,y2), ( 49 ) 

Xi 

where the asymptotic series for the hrst and second surface are given by xi = 
■51(1/1,1/2, \/h) and X2 = ■52(1/1,1/2, \/h)’ respectively. In addition, the series must satisfy 

dF dxi dF dxi dF 

- ^ - ^ ^ ^ 

oxi oyi 0x2 oy 2 oyi 

where i = 1 and 2, and F = Fq + pFi -|- /iF2 + •••, with Fq 7^ F^iy) but Fi, F2, ••• being 
functions of both x and y] see the text below Eq. (IT^ . 

Poincare assumed that there were certain values x^ for which dF/dxi are 
commensurable, which allows writing Xi = d)*(1/1,1/2) with F(<I>i, <I>2,1/1,1/2) = C 
satisfying Eq. fl 50 |) . To integrate this equation, Poincare considered Xi = x'^F (/ixj)^/^ -t- 
fvx] + ..., and determined the coefficients xf, with k = 1 , 2, .... The purpose of the second 
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Figure 7. Intersection of the asymptotic interf aces with the transver se section yi = 0 
as explained in the main text. Reproduced after iBarrow-GreenI 199^. Copyright 1994 
Springer. 


approximation was to evaluate an arbitrary number of coefficients of the above series, 
and to write the approximate equations for the asymptotic surfaces. Finally, in the 
third approximation he constructed the intersection of the asymptotic interfaces with 
the transverse section j/i = 0, and returned to his geometrical description, as depicted 

in Fig. [71 

In Fig. [71 the closed curve represented by the dotted and dashed line corresponds 
to the generating (unstable) periodic orbit. In addition, the curves AO'B' and A'O'B 
represent the asymptotic surfaces and their intersection with |/i = 0, and the dashed line 
represents the curve |/i = |/2 = 0. Poincare realized that the curves AO'B' and A'O'B 
could not be closed. He considered possible trajectories that would simultaneously 
belong to both si des of the asymp totic surface and called them doubly asymptotic or 
homoclinic orbits Andersonl . [1994 . 


Poincare identified primary, secondary, tertiary, and quaternary homoclinic orbits 
and realized that there could be an inhnite number of them leading to a homoclinic 
tangle, which appeared to be the hrst mathematical manifestation of chaos in the 
CR3BP; knowing how complicated the tangle could be, he did not attempt to draw 
it. Chaos has a sensitive dependence on initial conditions, and its existence in 
the CR3BP was confirmed by computational simulations of the CR3BP performed 
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Vela-Arevalo and Mardsen 


by ISmitli and Szebelielvi 1993 |. iKoleman et al.l 2012 | . and many others. Moreover, 


2004l | demonstrated that a method based on time-varying 


frequencies is a good diagnostic tool to distinguish chaotic trajectories from regular ones. 

Some of the modern co mputational me thods use Poincare sections to identify quasi- 
periodic and chaotic orbits. IPoincard jl892| studied trajectories close to a periodic orbit 
and investigated their return to a line drawn across the orbit (now known as Poincare 
section). The behavior of trajectories returning to the line forms the hrst return map, 
or Poincare’s map, and he used it to reduce the dimension of phase space of the three- 
body problem by one. He also showed that solutions could leave two equilibria and then 
return to them asymptotically, thus forming a heteroclinic orbit. The above concepts 
originally introduced by Poincare are iiow commonly used in modern theori es of chaos 


Thompson and Steward . Il986l . iHilbornl . 1 1994 iMnsielak and Musielakl . 12009 


7. The elliptic restricted three-body problem 
7.1. Basic equations and integrals of motion 

Similar to the CR3BP, the elliptic restricted three-body problem (ER3BP) requires that 
the two objects (called the primaries) have their masses significantly larger than the third 
one (Ml >> M 3 and M 2 >> M 3 ), and that the third mass has no gravitational influence 
on the primaries. However, the main difference between the CR3BP and ER3BP is that, 
in the latter, motions of the primaries are along elliptical orbits around their center of 
mass, which has important implications for the math ematical de s cripti o n of the s ystem 


as there is no t a Ja cobi’s constant in the ER3BP Szebeheld . 119671 . iMarchai 11990 


Palacian et ah. 2006 


In general, the third mass can move in 3D space, a case is called the 3D ER3BP. 
However, if motions of the third body are restricted to the same plane as the primaries, 
we have the planar ER3BP. Typically, the governing equations for the 3D ER3BP 


Szebehelv. 

1967. 

Marchal. 

1990 


In the following, we 

describe the governing equations of the 3D ER3BP, but then we make clear distinctions 
between the 3D and planar ER3BP when discussing some recently obtained results. 

Let M 2 move around Mi on an elliptic orbit with eccentricity e, true anomaly / 
and semimajor axis a = 1, and let y, rj, ^ be the rotating-pulsating coordinates. With / 
being an indepe ndent vari a ble, t he set of governing equations describing the 3D ER3BP 
can be written Szebehelv . Il967j | in the following form: 

d‘^X 


dP 

d'^T] 

w 


2^-n 

df~ ^ ’ 
dx 

+ ~ ’ 


(61) 


(52) 


, 21 
dp + df 




and 


( 53 ) 
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where w = ^2/(1 + ecosf) with ^2 = (x^ + + C^)/2 + (1 — /u)/ri + /x/r 2 + /i(l — /i)/2, 

= (x — + ? 7 ^ + and rf = (x — /i + 1)^ + For the planar ER3BP, it is 

necessary to take ^ = 0. _ 

We now follow Llibre and Pinol |l990 | and introdnce qi = —x + /i, ^2 = —Vi Qs = ^i 


Pi = —x' + Vi P '2 = —v' + P ps = where prime indicates d/df. Then Eqs fl5Tl) 
through fl53l) can be written as 


dqi OH 
df dpi 


dpi 

df 


OH 

dqi 


(54) 


where i = 1, 2 and 3, and H is the time dependent Hamiltonian H = Hq fxHi with 


-^0 — 2[iP^ + ^ 2 )'^ + (P2 — <?i)^ + pI + qI] 


1 


1 + e cos / 


and 


Hi = 


1 2 1 1 

TW- f ( - 

1 + e cos j \ ri r 2 


2 ( 9 ? + 92 + + ~ 


, (65) 


(56) 


Since H is an explicit function of /, it is not an integral of motion for the ER3BP; 
however, H is an integral of motion for the CR3BP for which e = 0. The planar ER3BP 
requires that qs = ps = 0. 

We may now apply the procedure of deriving the Jacobi integral for the CR3BP (see 
Section 6.3) to Eqs fl5T]) through ([53]), which are the ER3BP governing equations. The 
procedure requires multiplying Eqs fl5Tl) . fl5^ . and fl53l) by x^ v'i respectively, 

adding them and integrating. The result is 


df) \df 


+ |2i 

df 


= 2 


{u^dx + uf^df] + u^d^) , 


(57) 


which is equivalent to the Jacobi integral in the CR3BP; however, the difference is that 
the RHS of the above equation is not constant. Thus, the Jacobi integral for the ER3BP 
is not a constant of motion. 

There are othe r constants of motions for the ER3BP. As originally pointed out by 
ContopoulosI 1967( |. two integrals of motion for the planar ER3BP exist if the system is 
considered in rotating coordinates with their x -axis p a ssing through the primaries. The 
results were generalized to the 3D ER3BP by ISarrisI 1982| , who found three integrals 
of motion for a small eccentricity of the relative orbit of the primaries, and for a small 
distance of the third body from one of the primaries, and demonstrated that they 
depended periodically on time. 


1.2. Periodic solutions, libration points and Lie series solutions 

Periodic orbits show ing chains of collision orbits of the original Kepler problem 


Searches for such orbits in the 3 

3 ER3BP were performed by 

Gomez and Olle 

1991al. 

Gomez and Olle 1991b 

. Bertott 

1991 

. Bolotin and MacKav 

20 QQ 

. Bolotin 

2005 

and 
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Palacian and YanguasI 2006| . In Be rtotti’s work, a variational method designed to find 
such periodic orbits was developed. iBolotinI 2005| considered the planar ER3BP with a 


small mass ratio and eccentricity. Then he proved the existence of many second species 
periodic orbits. The fact that there are ejection-collision orbits in the 3D ER3BP was 


demonstrated by lLlibre and Pinoll 1990|. Moreover, the effectiveness of the phenomenon 


of the gravitational capt ure of the third body by the primaries in the 3D ER3BP with 
sma ll u was discus sed by iMako and SzenkovitsI [2004 . 


Sinclairl 197nl | identified and classified various families of periodic orbits in the 


commensurable planar ER3BP, where two bodies are close to commensurability in their 

id the central body. Moreover, 
discovered new families of stable 


Hadiidemetriou 

1992 

and 

HaehiehiDour et ah 

2003 


3:1 and 1:2 periodic orbits in the planar ER3BP. 

In addition to the existence of periodic orbits in the 3D and planar ER3BP, there 
are also the equilibrium point solutions or the libration points in both problems. The 
conditions that define the libration points are obtained from Eqs flHTll through fl53ll . 
namely, y" = rj" = C = x' = v' = C = ^ and dfl/dx = dfl/dr] = dfl/d^ = 0. Using 
these conditions, it is easy to show that there are five libration points in the ER3BP 
(similar to the CR3BP), a .nd that t hree o f them are collinear and two are triangular or 


equilateral iDanbvI. 


1964 iBennettl. Il965| . Motions around these libration points were 


investigated by ISzebehelyl jl967j | , who demonstrate d that the triangu lar points pulsate 
together with their own coordinate systems, and bv IChondhrvI 1977|, who showed that 
the coordinates of the libration points depend explicitly on time. 

who used the 
The 


A solution to the planar ER3BP was constructed by iDelva 


method of Lie series that is another form of the Taylor series [Schneiderl . 11979 


1984 


ER3BP in the rotating-pulsating coordinate system was considered and the Lie operator 
for the motion of the third body was d e rived as a function of coordinates, velocities, and 
true anomaly of the primaries Delv4 Il984| . An analytical form of the Lie series was 


obtained and used to compute orbits; an extension of this method to the 3D ER3BP 
should be straightforward. Numerical solutions to both the 3D and planar ER3BP 
were obtained and some specific examples are presented in the next section devoted to 
applications. 


7.3. Elliptic Hill problem 

We now briefly discuss the elliptic Hill problem and its potential applications to 
astrodynamics and extrasolar planetary systems. As stated in Sections 1 and 6.7, in the 


we then have the elliptic Hill problem Ichtiaroelou 

1980 

. A few families of periodic 

orbits in the elliptic Hill problem were discovered by ] 

chtiaroErlon 

1981 

, who found all of 
calculated a large 

them to be highly unstable. In more recent work. 

Vovatzis et ah 20121 


set of families of periodic orbits in the elliptic Hill problem, dete rmined their stability 


and applied the results to motions of a satellite around a planet. ISzenkovits and Makd 
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2008| established a more accurate criterion of the Hill stability and used it to investigate 


stability of exoplanets in the systems: 7 Cephei Ab, Gliese 86 Ab, HD 41004 Ab, and 
HD 41004 Bb. They concluded that the orbits of these exoplanets will remain stable, 
according to the Hill stability criterion, in the next few million years. 


8. Applications of the restricted three-body problems 

Among the possible solutions of the three-body problem, the simplihcation imposed 
in the restricted three-body problem has been the most practical, and therefore many 
different astronomical systems have been studied within this framework. We present 
applications of the restricted three-body problem with an emphasis on extrasolar planets 
and moons. Since orbital dynamics plays a role in determining the likelihood of a planet 
being classihed as habitable, residing in a region of space in which liquid water could 
persist on the surface of the planet, we provide additional discussion within that context. 


8.1. The Earth-Moon system with a spacecraft 


The advent of the space exploration program presented a host of challenging three-body 
problems, and it established a new area of study in astrodynamics. The initial problem 
was that of solving the motion of artihcial satellites, followed by the task of determining 
spacecraft orbital maneuvers required for missions to the Moon. 

The orbital dynamics of artihcial satellites and spacecraft is based on n-body theory, 
but it can be approximated as a three-body problem. Moreover, we distinguish between 
spacecraft trajectories and motion of artihcial satellites (orbiters, space telescopes, 
space stations, communications satellites), as the solution for each case have diherent 
requirements and outcomes. The following are some interesting examples: 

(i) For analysis of space missions, such a spacecraft moving between two planets, 
or from Earth to the Moon, we may apply the restricted three-body problem. However, 
space missions through the farther regions of the Solar System are rather complex and 
require that spacecraft trajectories be designed by decoupling an n-body system into 
several three-body systems. For example, by decoupling the Jovian moon n-body system 
into several three-body systems, it is possible to design a spacecraft orbit which follows a 
prescribed itinerary in its visit to Jupiter and its moons. In fact, the Jupiter-Ganymede- 
Europa-spacecraft s ystem is approximated as two coupled planar three-body systems 


Gomez et al.l . DOnij . Furthermore, for a spacecraft in transit from Earth to the outer 


planets, its flight may involve a gravity assist, which yields a four-body problem; once 
far away from the Earth, it reduces to a three-body problem. Further complications 
arise when spacecraft move in the vicinity of asteroids or other celestial objects with 
gravitational helds that are not sufficiently well-established. 

(ii) For artihcial Earth satellites, the problem of the orbit evolution is an instance 
of the four-body problem. The evolution of a satellites orbit with a high eccentricity, 
such as that of the Milniya communication satellites, is very different from that of a 
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low-Earth orbit of low eccentricity (e.g., Hubble Space Telescope). In the latter case the 
most important perturbations are due to the oblateness of the Earth and the air drag. 
Both effects are small for a moderate to high eccentricity orbit with a large semi-major 
axis. But the gravitational influence of the Sun and the Moon becomes predominant. 
The Kepler Space Telescope orbits the Sun, which avoids the gravitational perturbations 
and torques inherent in an Earth orbit. 

In general, the design of space missions requires analytical and numerical methods 
to solve planar circular restricted three-body problem (PCR3BP) and circular restricted 
three-body problem (CR3BP). For example, when decomposing an n-body system 
into three-body systems that are not co-planar, such as the Earth-Sun-spacecraft and 
Earth-Moon-spacecraft systems, we need to consider three-dimensional conhgurations 
and solve circular restricted three-body problems (CR3BP). The planar motion of a 
spacecraft relative to two other masses which orbit each other is very complex, with 
nonlinear effects due to the perturbation of the two bodies, for which no general 
analytical solution exists. Therefore, the orbital analysis for such a spacecraft is based 
on sophisticated numerical methods. 


8.2. The Sun-Earth-Moon system 


This particular solution of the three-bo dy problem has influenced the theory of gravity 


since its inception by iNewtonI 1687l |. In comparison to the previous application 


concerning satellites, this approach relies on the Sun-Earth system as the two dominant 
bodies with the Moon as a test particle. This was effectively the approach taken 
by Newton; however, he was not successful in determining stable solutions (~ 8% 
uncertainties). The problem is that the mass of the Moon is comparable to that of the 
Earth within two orders of magnitude; it is also relatively close to the Earth, therefore, 
inducing reaction forces and tides that are not accounted for in the restricted three-body 
problem. 

Others hav e made progress i n determining stability conditions for the Sun-Earth- 


Moon system. iHilll jl877l . Il878| developed a general stability criterion based upon 


the balance of gravitational forces on a test particle between the Earth and the 
Sun. This criterion has been rehned statistically through numerical simulation with 
estimates indicating a limit of ~ 3.2R h. where Rh denotes the Hill Radius given by 
Rh = a Glad man . 1993l |. Limitation s of this criterion in appli cations two 


extrasolar p l aneta ry systems were pointed ou t bv ICuntz and Yeager l2009l . Recently, 


Satyal et al.l |2ni3| used a method proposed bv iSzenkovits and Makdl [20081 1 that utilizes 
the Hill stability criterion; we discuss this approach in Section 8.6. 


8.3. The Sun-Jupiter system with an asteroid 

Beyond the Sun-Earth-Moon system, the next most studied restricted three-body 
problem is the Sun-Jupiter-asteroid system. Since masses of asteroids are very small 
when compared to either the Sun or Jupiter, their description hts well within the 
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framework of the restricted three-body problem. Specihc applications include dynamics 
of asteroids within the Asteroid Belt, and the origin of the Kirkwood Gaps related to a 
depletion of the asteroid population in some orbits due to re sonances with J upiter . 


2002 


Resonances of asteroids with Jupiter have been studied by iNesvornv et ah 
including possible three-body resonances in the Solar System. 

As discussed in Section 6.2, the collinear Lagrange points are unstable. Thus, no 
population of asteroids in the vicinity of these points is expected. However, there is a 
population of asteroids with stable orbits around the L 4 and L 5 Lagrange points. These 
populations are the Trojan asteroi ds, and their stability has been investigated usi ng the 


restricted three-body formulation Bentlerl. 120051 IValtonen and Karttnuenl . 12006 


8 .4- A single star with giant and terrestrial exoplanets 

Since the conhrma t ion o f the hrst exoplane t 51 P egasi around a solar-type star 


Mayor and Qnelozl . 119951 iMarcy and Butleii . Il995| . we now know almost 1800 


exoplanets, including more than 1000 exoplanets discovered by the NASA’s Kepler space 
telescope; this number of discovered planets is growing everyday 0. Prior to the launch 
of Kepler, most known exoplanets were Jupiter-analogs with respect to their sizes and 
masses. A very small subset of these systems are now known to host terrestrial planets. 
Thus, the restricted three-body problem has been used to investigate a planetary system 
with a single star, gas giant, a nd terrestrial exopla nets. 


In studies performed by iNoble et al.l 2002|, broad regions of orbital stability 


were found for terrestrial exoplanets located in the habitable zone (HZ) of 51 Pegasi 
and other planetary systems. More general stabi lity limits for extra-so l ar pla netary 


systems with two exoplanets were established by iBarnes and Greenberd 2006| . The 


HZ’s for main - seque nce stars of different spectral stars were originally established by 


Kasting et al.l 1993| . Prospects for the existence of a terrest r ial pl anet in stellar HZ 
with known giant exoplanets were discussed by Hones et ahl 2005| . In more recent 


work, iKopparapn and BarnesI 2nifll | identihed 4 extra-solar planetary systems that can 
support a terrestrial planet in their HZ’s. _ 


In the recent work of iGoldreich and Schlichtingi 2ni4| . resonances in the planar. 


circular, restricted three-body problem with possibilities of semimajor axis migration 
and eccentricity damping result from planet-disk interactions. The authors showed that 
the migration and damping results in an equilibrium eccentricity, and that librations 
around this equilibrium are overstable. Applications of these results to exoplanets 
discovered by the Kepler satellite demonstrated that most of the observed planets are 
overstable, which leads to passage through resonance. A more realistic case of the 
planar, three-body problem is also considered with two (inner and outer) exoplanets 
being in resonance, and the departure from exact resonance is discussed. 

The most prolihc detection method before Kepler was the radial velocity method, 
which relies on shifts of lines within the spectra of the host star. These shifts demonstrate 


^ http://exoplanet.eu/catalog/ 
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the reflex motion of the star due to the gravitational force exerted by the unseen 
companion. Terrestrial exoplanets would induce shifts in the radial velocity of the 
host star on the order of centimeters per second, which is just beyond today’s limits 
for this detection technique. In terms of exoplanet detections, the photometric transit 
method has uncovered through the it'ep/er sample, that giant and terrestrial planets can 
exist wi thin the same system with a larger fraction of planets being Neptune-sized or 

In such systems 


smaller Howard et al. 


2012 


Fressin et al. 


2013 


Santerne et al. 


2013 


the restricted three-body problem can be used to constrain the likelihood of additional 
planets. 


8.5. A single star with a giant exoplanet and exomoon 

Another set of sy stems has been proposed for the Jupiter-like planet s discovered within 


a host star’s HZ Kasting et al 


1993 


Kopparapu and Barnesl. 


2010 


Giant exoplanets 

within this region present a major challenge to the search for terrestrial exoplanets, 
which could support life. It has been shown that giant exoplanets near the HZ can 
induce perturbations upon the putative terrestrial exoplanet, thereby making the orbit 
of the smaller body unstable. In order to keep the search for habitability meaningful, 
exomoons have been considered as possible abodes for life, which orbit giant exoplanets 
within the HZ. Studies of such systems are based on the restricted three-body problem. 

le habitable zone to host 
The formation of 


The plausibility of a Jupiter-like exopl anet located in t 


an Earth-mass moon has been considered Williams et 


al. 


the capture of a terrestrial planet still remains possible Kipping . 120091 . iKipping et al. 


1997 


such a system may experience considerable difficulties Canup and Wardl. 


2006 


but 


2012, 2013b 


a 


Other studies have used a parameter survey approach to determine 
the stability of pos sible exomoons; see applications to the system HD 23079 system 


Cuntz et al.l . l2013| . Studies of such systems are based on the formulation of the 


restricted three-body problem. 

8 .6. Binary stellar systems with a giant or terrestrial exoplanet 

One of the most extreme cases of the restricted three-body problem involves the study 
of exoplanets orbiting either one or both components of a binary star system, in 
which the more massive and less massive stars are called the primary and secondary, 
respectively. The categorization of these systems has been referred to as satellite-like 
(S-type) or planetary-like (P-type), where the former contai ns a planet that orbits 


a single star and the latter orbits both stars Dvorak . 1982 1 ; the extent of HZ’s in 


such systems have been recentl y estab li shed [Eggl et al 


Kaltenegger and Haghighipourl . 120131. iHaghighipour and Kalteneggerl . 


2012 


Kane and Hinkel. 


20131 . 


2013 


Cuntz], 


2014| . Using this convention, many different systems have been studied based on 
the restricted three-body problem, with various methods including topology and chaos 
indicators. Furthermore, planetary formation in these systems have been studied 
even before one P-type system (Kepler 16) with the exoplanet had been rigorously 
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Figure 8. Determination of stability of S-type planets using the method of Lyapunov 
exponents. Regions of stability (blue) are indicated as well as regions of chaos 
(red). Areas with a white background denote unstable regions on 0.1 Myr timescales. 


Reproduced after lOuarles et al.l [2011j . Copyright 2011 Astronomy & Astrophysics. 


confirmed by 

Dovle et ah 

2911 

Pilat-Lohinger and Funk. 

2 niol|, 


, As this field encompasses a separate research area 
we now discuss only selecte d rece nt approache s. _ 


Using a topological method, iMusielak et ah 


2005 


and Cuntz et ah 2007 


determined stability criteria for the special case of binaries hosting exoplanets, which 
relied upon the Jacobi constant formalism in contrast to previous statistical approaches 


Holman and Wiegerti llOOOL iDvorakl. Il984| . This approach allowed for the exploration 


of marginally stable systems, which would have prematurely been characterized as stable 
using a parameter space consisting of the mass ratio p, = M 2 /(Mi + M 2 ) and the initial 
distance ratio po = R/D, where M 2 < Mi, R is the distance from the planet to Mi, and 
D is the initial separat i on of t he stellar components. 


Eberle and Cnntzl |2nina| furthered this exploration of d etermining stab ility of 


planets in stellar binary systems using Hamilton’s Hodograph Hamiltonl. Il847l |. This 


method demonstrated how the effective eccent ricity of an exoplanet wou ld change due 
to interactions with both stellar components. lEberle and Cnntzl 2ninbl | explored the 
same parameter space and produced a ma p of general stability with contours of the 
effective eccentricity. iQnarles et al.l 2nil| use d the parameter space along with the 


maximum Lyapunov exponent Twapnnovl . Il907l |. and found similar contours of stability 
with structures relating to orbital resonances and islands of instability. 

Figure | 8 ] s hows the parameter space with stability contours within a linear regime. 


Quarles et al.l 201 1 | showed that a 3:1 resonance can cause chaos within the circular 


restricted three body problem for a wide range of mass ratios. Although these results 
have been mainly theoretical, they have been applied to the limited number of exoplanets 
found in stellar binary systems. One such system, z/ Oct, has a controversial discovery 
as the proposed distance ratio of the exoplanet lies within a region of instability (see Fig. 
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0.36 0.38 


0.36 0.38 



Figure 9. Determination of stability of the planet in ly Octantis assuming a retrograde 
orbit ing planet initially p laced at apastron (left) and periastron (right). Reproduced 


after lOuarles et al.l 2012a| . Copyright 2012 Monthly Notices of the Royal Astronomical 
Society. 


9). This posed a problem for the conhrmation of the exoplanet. However, if one considers 
an exoplanet in retrograde orbit relative to the orbiting binary, then enlarged regions of 


2010 b. 

Quarles et ah. 

2012 a. 

Gozdziewski et ah. 

2013 


The rest ricted three-body p roblem has also been applied to P-type or circnmbinary 


exoplanets. iDovle et al.l 201 1| annonnced the hrst circnmbinary exoplanet in the 


Kepler-16 system. This has been an important discovery for the theory behind the 
restricted three-body problem to be tested. Soon after the discovery was annonnced. 


Quarles et ahl 2012bl | investigated cases when an Earth-mass exoplanet could be injected 


into the system; the studies of orbital stability of such an exoplanet could help in 
discovering habitable Earth-analogs in the Kepler-16 system. 

Since injected terrestrial planets in the Kepler-16 system constitutes a restricted 


four-body system, iQuarles et ahl |2012bl | reduced it to a restricted three-body problem 


in order to determine the possibility of the existence of Trojan moons in this system. 
With the distance between the stellar hosts (denoted as A and B in Fig. 10), and 
the giant (Saturnian) exoplanet (denoted as b in Fig. 10) in the Kepler-16 system 
being relatively large. Earth-mass objects could remain in stable orbits around the 
approximated equilateral (L4 and L5) points in many different types of orbits. Note 
that the giant planet is located at the border between the HZ and the so-called extended 
habitable zone (EHZ), wh ich requires rnore ext reme exoplanetary atmosphere with a 


signihcant back-warming Quarles et al.l. l2012b| . Figure [TO] shows one of the orbits 


possible in this sys tem; an exop lanet in this regime could be habitable, if it has a strong 
greenhouse effect CnntTi . l2014l |. However, either the formation or capture of such an 


exoplanet would be difhcult to achieve with in the environment of the circnmbinary disk 


Meschiari. 

2012 , 

Paardekooper et ah. 

2012 

























































The three-body problem 


45 



Figure 10. Orbits of possible Trojan moons in Kepler-16 using a three-body 


approx imation - see the main text for explanation. Reproduced after iQuarles et al 


2012bl| . Copyright 2012 The Astrophysical Journal. 


9. The relativistic three-body problem 


9.1. The Newton and Einstein theories of gravity 

The three-body problem discussed in the previous sections is based on the Newtonian 
theory of gravity in which the bodies are assumed to be spherical and can be considered 
as point particles of given masses. Since Newton’s description of gravity is only a weak 
held approximation of Einstein’s General Theory of Relativity (GTR), which is currently 
our best theory of gravity, now we want to discuss implications of using GTR on the 
three-body problem. GTR replaces the hat space-time of Special Theory of Relativity 
by a 4D smooth, continuous, pseudo-Riemannian manifold endowed with the metric 
ds^ = g^^{x)dx^dx'^, where g^y is the metric tensor to be determined by Einstein’s held 
equations, u and u a re 0, 1, 2 and 3, and the usual summation conventions are employed 


Hobson et al.l. l2006j . According to GTR, the 4D space-time is curved by the presence 


of bodies and the resulting curvature is identihed with gravity, which clearly shows that 
Newton’s and Einstein’s theories of gravity are signihcantly diherent. Nevertheless, 
Einstein developed his GTR in such a way that in the limit of weak gravitational helds, 
Newton’s description of gravity remains valid. 

To account for relativistic ehects on motions of astronomical bodies, relativistic 
celestial mechanics was formulated and its methods appli ed to many inte r esting 


astronomical problem s, including the three-body problem Brnmbergl . Il972[ Il991 


Kopeikin et al 


2011 


Since Einstein’s held equations of GTR can only be solved for the 
simple one-body problem, and even the two-body problem cannot be solved rigorously 
in GTR, the equations of motion of the relativistic three bodies must be obtained 
with some approximations. Typically, it is the so-called post-Newtonian approximation, 
which is a weak held and slow motion approximation, that is used to derive the post- 
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Newtonian equations of motion of the three spherically symmetric bodies; the resulting 
equations depend not only on masses but also on physical parameters relevant to the 
interna l structure of th e bodies, such as elastic energy, mo ment of inertia, among 


others Brumberd . Il99l| . According to iKopeikin et ahl 201 1| and references therein. 


the presence of these parameters on the motions of the three bodies may actually 
be irrelevant. Thus, some authors prefer to use the Einstein-Infeld-Hoffmann (EIH) 
equation of motion, which is va lid for point-like masses with no internal structure 
Einstein et al. . 1938 . Infel . [l9^ . 


9.2. The general relativistic three-body problem 

In the general relativistic three-body problem, there is no restriction on masses, 
except that they are hnite. The po st-Newtonian equations of motion for this problem 
were obta i ned b y iBrumberd 1972| and applied t o different astronomic al settings by 
Brumbergi 199l[ | as well as to the Solar System by iKopeikin et al.l 2nil| . Hamiltonian 


and standard formalisms have also been forn iulated for the cons i derat i on of genera l 
relati vistic effects within an n-body framework Saha and Tremainel . 119921 1 19941. iKidderl . 


1995 


The existence of the post-Newtonian collinear solution, which is a relativistic 
extension of the original Euler colli near solution (see Section 3.1 and Fig. 1), was 
shown by lYamada and Asadal 2010al | who used the EIH equation of motion. They also 
calculated the relativistic corrections to the spatial separation between the three bodies 
required for the soluti on to be valid. The uniqueii ess of the relativistic collinear solution 
was demonstrated by lYamada and Asadd j2ninb . 


Moreover, a triangular solution that is a relativistic version of the Lagrange solution 
with an equilateral tri angle ipi n ing th e three bodies (s ee Section 3 . 1 and Fig. 2) was 


originally obtained by iKrefetzl 1967j. More recently, llchita et al.l 2nil| studied the 


post-Newtonian effects on this solution and showed that the solution satishes the post- 
Newtonian equat ions of motion only when a ll three masses are equal. Thereafter, it was 
demonstrated bv lYamada and Asadal 2012l | that the relativistic corrections to each side 
of the triangle must be added and, as a result, the relativistic triangle will not always 
be equilateral. 

Another interesting result is that the 8-type periodic solution for the general (non- 
relativistic) three-body problem with equal masses (see Section 3.2 and Fig. [3|i) does also 
exist in the general relativistic three-body problem. Its hrst a nd sec ond post-Newtonian 


orders were investigated numerically by llmai et al.l 2nn7l | and iLousto and Nakano 


2008 


respectively, who demonstrated that the relativistic effects are responsible for 
only small changes in the shape of this 8-type solution. 

It is also likely that re l ativis tic versions of the 13 new periodic orbits found by 


Suvakov and Dmitrasinovid 2013| in the planar general (non-relativistic) three-body 


problem (see Section 3.2 and Fig. [3]) will be discovered and that relativistic corrections 
to the results presented in Fig. |3]will be found. However, at the present time this idea 
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9.3. The restricted relativistic three-body problem 


In the restricted relativistic three-body problem (RR3BP), the two primaries have 
dominant masses and move around their center of mass; however, the third mass is very 
small and its gravitational influence on the primaries is negligible. The post-Newtonian 
eq uations des c ribing motions of the third body in the RR3BP were originally obtained 
bv iBrumberd 1972l |. who used a synodic (rotating) frame of reference with the origin 
at the center of mass and the primaries fixed. 

L et us consider a planar RR3 BP, introduce the synodic coordinates v and C, and 

and write the post-Newtonian equations of motion 

dU 


follow iDouskos and PerdiosI 2002 
in the following form 


i) — 2nC, = 


dv 


d (dU 
dt \ di) 


and 


- dU d fdU 

c 2 nv = —^ 

dC dt\ dC 


(58) 


(59) 


where ‘dot’ represents the derivative with respect to time t, n = 1 —3[1— /r(l —/i)/3]/2c^, 
with c being the speed of light. In addition, U = Uc+Ur is the RR3BP potential function 
composed of the classical potential 74 = r‘^/2 -\- (1 — p)/ri -\- pi/r 2 , with r = \A4”+T^, 
ri = \/{v -\- pbY -\- and r2 = ■^/(u -|- /i — 1)^ -|- The relativistic corrections 74 are 
given by 


74 = —6r^c^ 1 — -/i(l — /i) 

o 


c" 


6c^ 


1 + n 

ri r2 


— 2c^p,(l — fi) 


n 


(v +cr +(c - 

(1-3/i 


+ 0^ + (C - V? 

1- fi 


2 c^ 


n 


r 2 




ri r2 / 


80 + CM3 + 


(60) 


The libration points for the RR 3BP are obtained by s ubstituting i; = C = 'b = C = 0 
in the above equations of motion Contopoulos . 1976 . Bhatnagar and Hallan . Il998| . 
These libration points are the relativistic counterparts of the Lagrange collinear Li, L 2 
and L3, and triangular L4 and point s . Lin ear st ability of the relativi s tic co llinear 


points was investigated by iRagos et al.l j2nnn| and iDouskos and PerdiosI (2002|, who 


showed that all these points were unstable, which is consistent with the results obtained 
for the non-relativisti c collinear points. 

linear stability of the relativistic 


In the work by iBhatnaear and Hallan 


1998 


triangular L4 and L5 points was studied and it was demonstrated that these points were 
unstable for the whole range 0 < p < 0.5, despite the well-known fact that the non- 
relativistic L4 and L5 are stable for p < po, where po = 0.038521 is the Routh critical 
mass ratio (see Section 6.2). The problem was later revisited bv IDouskos and Perdios 
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2002 l | and lAhmed et al.l |2006| , who found that the relativistic triangular points are 
linearly stable in the r ange of mass rations 0 < u < Ur, wher e Ur = fto — 17\/^/486c^ 


Donskos and Perdiosl . |2002|, and j^r = 0.03840 Ahmed et ahl . 12006 


Amon g different possible ap plications of the RR3BP as discussed by Brumber^ 


1991l | and iKopeikin et al.l |2011l 


of Mercury’s perihelion made by iMandl and DyorakI [1984 


let us also mention the calculation of the advance 

the computa tion of 


chaot ic and non-chaotic trajectories in the Earth-Moon orbital system by IWanex 
the relativistic corrections to the Sun-Jupiter libration points computed 

the analysis of stability of c ircular orbits in the 


2003 


by lYamada and Asada 


2010 a 


Schwarzschild-de Sitter space-time performed by lPalit et al.l j2009| . and the investigation 
of the GTR effects in a coplan ar, non-resonant planetary systems made by 


Migaszewski and Gozdziewskil 2009 


10. Summary and concluding remarks 

We haye reyiewed the three-body problem in which three spherical (or point) masses 
interact with each other only through gravitational interactions described by Newton’s 
theory of gravity, and no restrictions are imposed on the initial positions and velocities. 
We began with a historical overview of the problem, where a special emphasis was 
given to the difficulty of Ending closed form solutions to the three-body problem 
due to its unpredictable behavior. We gave detailed descriptions of the general and 
restricted (circular and elliptic) three-body problems, and described difierent analytical 
and numerical methods used to find solutions, perform stability analyses, and search 
for periodic orbits and resonances. We also discussed some interesting applications to 
astronomical systems and spaceflights. In the previous section, we also presented the 
general and restricted relativistic three-body problem and discussed its astronomical 
applications. 

The three-body problem described in this paper may be considered as classical 
(or standard) because the bodies were assumed to be spherical objects or points of 
giyen masses as described by either Newton’s theory of gravity (Sections 1-8), or by 
Einstein’s theory of gravity (GTR) in its post-Newtonian approximation (Section 9). 
However, there are studies of the ef fects of oblateness as well as Goriolis and centrifugal 
forces on the three-body problem jAbouelmagd et ah . 2013], and the existence of the 
libration points i n the restricted three-b ody problem with one or both primaries being 


Sharma et af 


2001 


oblate spheroi ds Markellos et al.l . Il996| . or both primaries being triaxial rigid bodies 


Moreoyer, the existence of periodic orbits in the restri cted 


three-body problem with one primary being an oblate body Mittal et al.l . Il996 


or 


with all objects being oblate spher oids and the primaries being radiation sources 
jAbouelmagd and El-Shaboury . 2012], have been investigated in the elliptic restricted 


three-body problem. 

Note that we 

cite only selected papers in which references to previously published papers on these 


Kumar and Ishwar. 

2011 , 

Singh and Umar. 

2012 
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topics are given. 

Also of interest is the so-called Chermnykh problem in which a point mass 
the gravitational field produced by an uniformly rotating dumb-bell 


moves m 


Gozdziewski and Macieiewskil . Il998l. iGozdziewskil . |2003| . This problem generalizes two 


classical problems of celestial mechanics, namely, the two-hxed center problem and the 
restricted three-body problem, with the former applied to the semi - classic al quantization 


of the molecular ion of hydrogen jGozdziewski and Macieiewski . 19981; however, any 


applications of the three-body problem outside of celestial mechanics and dynamical 
astronomy are outside of the scope of this paper. 

We would like to point out that our reference list is only a small subset of all 
papers and books published on the three-body problem. The choice of references was 
determined by what we considered to be the most relevant to the topics covered in this 
paper; however, we do respect the fact that other authors may have different opinions in 
this matter. Moreover, our reference list contains some papers with astronomical data 
where the results directly relate the applications of the general and restricted three-body 
problem to the Solar System, and newly discovered extrasolar planetary systems. 

We hope that we presented a balanced view of the three-body problem, and that 
our choice of specihc topics covered in this paper represents a good recognition of the 
previous and current research in the held. We also do sincerely hope that this review 
may serve as a guide to the major previous and current achievements in the held, and 
as an inspiration to scientists and students for opening new frontiers of research in this 
truly remarkable and very interesting problem. 
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